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The  results  of  the  research  effort  include  four  parts: 

1.  a)  The  Effect  of  Cementation  on  the  Elastic  Properties  of  Gran¬ 
ular  Material  and  b)  Seismic  Velocities  in  Compacting  Sediments 

2.  a)  The  Mechanics  of  Hollow  Grains  and  b)  Pillow  Basalts  Com¬ 
paction  due  to  Grain  Collapse 

3.  Modeling  the  Granular  Structure  of  Rocks 

4.  Anisotropic  Poroelasticity  and  Biot’s  Parameters 

PART  1  A  simple  analytical  model  has  been  developed  to  describe  the 
mechanics  of  cemented  granular  material,  for  which  a  solution  has  been  ob¬ 
tained  to  model  the  normal  and  shear  deformation  of  a  cement  layer  with 
straight  boundaries.  The  problem  of  the  normal  compression  of  two  con¬ 
tacting  circular  elastic  grains  and  an  elastic  cement  layer  between  them  has 
been  reduced  to  a  linear  integral  equation.  This  problem  has  been  solved 
both  for  two-dimensional  (cylindrical  grains)  and  three-dimensional  (spher¬ 
ical  grains)  cases.  The  elastic  modulus  of  such  a  combination  is  strongly 
increased  by  cementation  compared  to  the  case  of  a  classic  Hertzian  con¬ 
tact.  The  stiffness  of  the  cemented  system  increases  with  the  length  of  the 
cement  layer  and  with  the  relative  stiffness  of  the  cement.  It  significantly 
exceeds  the  stiffness  of  the  Hertzian  system  at  low  confining  pressure.  A 
small  increcise  in  the  amount  of  cementation  results  in  significant  growth 
of  the  contact  zone  between  grains  and,  thus,  dramatically  increases  the 
stiffness  of  the  system. 

The  theory  of  cementation  is  used  in  a  theoretical  model  for  the  compaction 
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of  sediments  formed  by  quartz  grains  separated  at  their  contacts  by  a  ce¬ 
menting  component.  Compaction  occurs  due  to  lithostatic  pressure  that 
results  in  increasing  contact  stresses  with  burial.  As  the  stresses  at  the  con¬ 
tacts  become  larger,  the  cementing  material  yields,  the  separation  between 
the  grains  decreases,  and  the  radius  of  the  contact  becomes  larger.  This 
changing  geometry  gives  rise  to  an  increase  of  P-  and  S-  velocity  with  con¬ 
fining  pressure,  and  therefore  with  depth.  Seisnoic  velocities  are  evaluated 
using  standard  formulas  relating  the  elastic  properties  of  granular  material  to 
the  contact  stiffness  between  two  grains  and  the  average  number  of  contacts 
per  grain.  Velocities  increase  rapidly  during  the  initial  stage  of  compaction. 
As  compaction  continues,  the  velocities  gradually  approach  constant  values 
that  correspond  to  direct  contact  of  the  grains.  An  interesting  effect  of 
smooth  peaks  in  P-  and  S-velocity  before  they  reach  their  constant  values 
is  predicted  if  cementation  is  very  soft  compared  to  the  grain  material  (e.g. 
quartz  grains  in  silt).  This  is  due  to  the  fact  that  the  normal  stiffness  of 
two  spherical  grains  separated  by  a  thin  layer  of  a  soft  cement  can  be  higher 
than  that  of  directly  contacting  grains. 

PART  2  In  order  to  model  the  mechanical  behavior  of  pillow  basalts, 
we  have  developed  a  theory  for  the  two-dimensional  deformation  of  hol¬ 
low  grains  under  external  loading.  We  model  the  grains  as  cylindrical  shells 
of  a  closed  cross-section.  This  theory  allows  us  to  describe  the  deformation 
of  arbitrarily-shaped  grains  with  walls  of  varying  thickness.  Calculating  the 
normal  and  shear  stiffness  of  hollow  grains,  one  will  be  able  to  determine 
effective  elastic  properties  of  their  ^gregate  using  the  expressions  devel¬ 
oped  for  a  random  packing  of  spheres.  Not  only  can  this  theory  be  directly 
applied  to  estimating  the  acoustical  and  mechanical  properties  of  basaltic 
pillows,  it  is  also  appropriate  for  studying  the  behavior  of  pelagic  sediments 
when  these  are  composed  of  hollow  grains  (for  example,  diatomites  or  other 
microfossils). 

Assuming  that  the  compaction  of  these  materials  is  by  grain  collapse  due 
to  an  increase  in  confining  pressure  that  exceeds  their  strength,  we  next 
derived  simple  criteria  for  the  collapse  of  2-D  thick- walled,  and  thin- walled 
circular  grains.  These  criteria  can  be  used  to  estimate  the  depth  of  transition 
between  low-velocity  (high  porosity)  and  high  velocity  (low  porosity)  zones. 

PART  3  The  interaction  of  rock  grains  with  friction  and  slippage  is  modeled 
using  an  analytical  solution  for  plain  elastic  deformation  of  cylinders.  The 
model  explains  the  difference  between  “static”  and  “dynamic”  elastic  moduli 
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of  naturaj  rock.  The  numerical  experiments  performed  on  a  four-grain  unit 
predict  the  dependence  of  ‘‘static”  elastic  moduli  on  confining  pressure,  on 
the  friction  coefficient  between  grains,  and  on  their  initial  positions.  The 
“dynamic”  elastic  moduli  depend  only  on  the  relative  initial  positions  of  the 
grains  and  on  the  confining  pressure.  Four-grain  units  with  varying  initial 
configurations  and  friction  coefficients  can  be  used  for  modeling  granular 
material.  The  “static”  and  “dynamic”  moduli  of  this  material  will  depend 
on  the  statistical  distribution  of  the  above-mentioned  properties. 

PART  4  Prediction  of  wave  propagation  in  a  submarine  environment  re¬ 
quires  modeling  the  acoustic  response  of  ocean  bottom  sediments  which 
generally  consists  of  porous  granular  materials  partially  or  wholly  saturated 
with  water.  The  effect  of  anisotropy  has  to  be  incorporated  into  the  model 
in  order  to  simulate  more  realistic  responses.  Following  Biot’s  theory  we 
present  a  formulation  for  seismic  wave  velocities  in  a  general  anisotropic 
poroelastic  medium.  We  also  identify  the  anisotropic  parameters  that  need 
to  be  evaluated  or  measured  for  the  purpose  of  estimating  the  velocities. 

The  motivation  for  the  work  summarized  above  and  appended  below  has 
been  to  understand  the  physical  properties  of  materials  in  the  upper  few 
hundred  meters  of  the  shallow  ocean  floor  (both  sediments  and  basalts). 
Particularly,  we  wanted  to  study  the  relationship  between  physical  proper¬ 
ties  (e.g.  elastic-  wave  velocity  and  attenuation,  or  strength  under  load)  and 
the  state  of  the  material  (micromechanics,  porosity,  composition  and  mor¬ 
phology).  The  approach  we  chose  was  to  construct  simple  physical  models, 
isolating  these  parameters  to  individually  determine  their  effects  on  physical 
properties.  The  parameters  we  chose  to  concentrate  on  were  (1)  the  geom¬ 
etry  of  grains  and  (2)  the  characteristics  of  the  material  at  their  contacts, 
which  must  transmit  vibrations  and  static  stresses  within  the  material.  In 
general,  it  is  these  properties  which  control  both  the  initial  properties  of  the 
materials  as  they  are  deposited,  and  the  change  of  those  properties  during 
compaction,  lithification,  and  diagenesis.  Five  of  the  six  papers  we  append 
deal  specifically  with  the  development  of  theories  to  describe  the  effects  of 
cement  (Papers  la  and  lb)  and  grain  geometry  (Papers  2a,  2b,  and 
3).  The  last  paper  (Paper  4)  is  a  summary  of  Biot  theory,  which  was 
undertaken  to  provide  a  framework  for  future  work  in  relating  our  microme¬ 
chanical  models  to  a  generally  accepted  macroscopic  model.  This  summary, 
although  not  at  the  level  of  application,  is  important  in  defining  directions 
for  theoretical  work  to  be  undertaken  in  the  future. 


3 


The  application  of  the  theories  we  developed  in  the  past  year  to  the  prop¬ 
erties  of  sediments  are  (1)  to  understand  the  increase  of  velocity  with  depth 
of  sediments  due  to  compaction  and  prior  to  diagenesis,  and  (2)  to  investi¬ 
gate  whether  mechanical  effects  may  play  a  role  in  the  abrupt  increase  in 
velocity  often  observed  at  so-called  “diagenetic  fronts”.  After  developing  a 
model  for  the  elastic  properties  of  grains  with  cement  (Paper  la),  we  used 
this  model  to  calculate  the  velocities  of  wave  propagation  in  sediments  with 
three  principal  components:  grains,  cement,  and  fluid-filled  pore  space.  The 
state  of  the  sediments  is  controlled  by  the  average  number  of  grain  contacts 
and  the  properties  of  those  contacts.  Two  models  were  developed  to  study 
the  change  of  these  properties  with  depth.  The  first  was  a  simple  frictional 
model  with  which  predictions  could  be  made  regarding  the  behavior  of  gran¬ 
ular  materials  under  load  (Paper  3),  and  their  recovery  upon  the  removal 
of  load.  The  second  model  (Paper  lb)  allows  the  construction  of  an  initial 
material  (the  sediment  as  it  is  deposited)  and  the  prediction  of  its  behavior 
under  steadily  increasing  load  (as  it  is  buried  by  further  sediments).  These 
theories  are  complementary.  Both  also  allow  us  to  investigate  the  capability 
at  any  depth  of  such  materials  to  carry  static  loads  (for  example,  from  struc¬ 
tures  on  the  seafloor).  The  latter  model  predicts,  once  the  material  begins  to 
deform,  a  smooth  variation  in  velocity  with  depth  similar  to  that  presumed 
by  phenomenological  models  for  compaction.  However,  the  compaction  rate 
and  hence  the  rate  of  change  of  velocity  is  determined  by  the  strength  of  the 
material  between  the  grains.  As  this  strength  can  be  measured,  the  model 
can  be  tested  quantitatively. 

One  interesting  result  of  these  models  is  that  under  certain  circumstances, 
materials  with  a  small  residual  cement  layer  are  stiffer  than  materials  for 
which  the  grains  are  in  direct  contact  at  a  point.  Although  this  must  be  in¬ 
vestigated  more  carefully,  it  does  offer  an  intriguing  mechanical  explanation 
for  the  occasional  observation  of  a  high-velocity  “layer”  at  a  diagenetic  front. 
A  consequence  of  the  frictional  model  is  that  the  rate  of  change  of  veloc¬ 
ity  is  quite  modest  until  the  frictional  strength  of  the  contacts  is  exceeded, 
whereupon  a  rapid  compaction  and  consequent  rapid  velocity  increase  oc¬ 
curs.  Again,  this  mimics  the  characteristics  of  a  diagenetic  front. 

Obviously,  neither  of  these  overly  simple  models  provides  a  complete  answer. 
However,  each  provides  useful,  testable  insights  into  processes  that  must 
occur  in  compacting  sediments.  And  the  rates  of  change  of  velocities  in 
each  case  are  controlled  by  the  strengths  of  the  contacts  between  the  grains. 


whether  such  contacts  are  cemented  or  are  free  to  slip. 

Attempts  to  adapt  models  of  granular  materials  to  the  behavior  of  basalts 
have  previously  met  with  some  resistance  from  members  of  the  scientific 
community.  For  example,  volcanologists  and  others  feel  that  the  underlying 
physics  is  entirely  wrong.  This  is  largely  because  of  the  conceptual  hurdle 
that  must  be  overcome,  because  we  must  adapt  the  scale  so  that  pillows  are 
the  “grains**,  which  are  several  orders  of  magnitude  larger  than  grains  in  sed¬ 
iments.  Once  these  hurdles  are  overcome,  however,  the  analogies  developed 
through  these  exercises,  which  again  have  been  found  to  be  in  quantitative 
agreement  with  data,  can  be  quite  valuable. 

With  this  in  mind,  we  developed  a  model  for  the  behavior  of  a  cylindrical 
tube  with  arbitrary  cross-section  under  load  (Paper  2a).  Such  a  model 
was  selected  because  of  the  common  observation  that  pillow  lavas  often 
form  as  tubes  through -which  material  in  passing  leaves  a  frozen  sheU  and 
a  partially  evacuated  interior.  These  tubes  are  often  cut  by  radial  cracks 
generated  by  cooling  stresses.  Once  the  characteristics  of  these  tubes  have 
been  determined,  construction  of  an  aggregate  is  straightforward.  Again,  the 
aggregate  has  several  useful  properties.  First,  it  is  substantially  “weaker” 
than  the  basalt  itself,  leading  to  very  low  elastic-wave  velocities.  Second, 
it  is  capable  of  rapid  velocity  increase  (and  porosity  reduction)  at  a  depth 
corresponding  to  the  collapse  pressure  of  the  tubes  (Paper  3b.  One  of  the 
principle  results  of  this  work  is  that  there  is  no  way  to  predict  the  relation¬ 
ship  between  velocity  and  porosity  for  such  materials,  particularly  inasmuch 
as  the  properties  of  the  tubes  are  controlled  by  the  position  and  number  of 
radial  cracks.  This  model  complements  our  earlier  work  describing  the  ef¬ 
fects  of  grain  contact  stiffnesses  in  pillows  on  the  properties  of  the  aggregate. 
Thus  the  grain  contact  stiffness  results  can  be  incorporated  quite  naturally 
into  a  more  complete  model  for  pillows  using  as  a  model  for  the  grains  the 
thin-walled  hollow  cylinder  developed  under  this  contract. 

The  attached  summaries  and  preprints  deal  primarily  with  the  first  phase  of 
this  work  -  that  is,  development  of  the  models.  The  next  step  -  application 
of  these  models  to  sediments  and  basalts,  and  comparison  of  the  predictions 
to  real  data,  is  now  under  way.  This  requires  quantitative  measurements  of 
the  properties  of  the  components  of  the  aggregates.  Specifically,  we  need  to 
know  the  moduli  of  the  matrix  material  and  of  the  cement,  in  the  case  of 
the  elastic  properties  measurements,  and  the  yield  strength  of  the  cement 
in  the  case  of  the  plastic  model  for  compaction.  In  studying  the  collapse  of 


basalt  pillow  tubes,  we  need  to  know  the  strength  of  the  matrix  material 
(basalt).  These  data  are  in  some  cases  available  in  the  literature.  In  the 
case  of  basalt  strength,  however,  there  is  relatively  little  data,  and  these 
measurements  are  now  being  made  by  members  of  our  group. 

Over  the  next  several  months  these  results  will  be  incorporated  into  the 
models,  and  the  predictions  will  be  tested  quantitatively  against  in  situ 
measurements  of  velocity.  Two  papers  (one  invited)  will  be  presented  at  the 
American  Geophysical  Union  Fall  1991  Meeting: 

Dvorkin,  J.  and  D.  Moos  (1991),  Seismic  velocities  in  compacting  sediments 
(abstract),  EOS,  Trans.  AGU  72:  Fall  Meeting  1991  Program  and  Abstracts. 
Moos,  D.  and  J.  Dvorkin  (1991),  Modeling  the  seismic  properties  of  shallow 
oceanic  Crust:  Effects  of  extrusive  morphology  on  seismic  velocity  (Invited 
abstract),  EOS,  Trans.  AGU  72  Fall  Meeting  1991  Program  and  Abstracts. 

Other  publications  acknowledging  support  from  this  contract  include: 

Dvorkin,  J.,  Mavko,  G.  and  Nur,  A.  (1991),  The  effect  of  cementation  on  the 
elastic  properties  of  granular  materials,  in  press  (Mechanics  of  Materials). 

The  remaining  papers  attached  to  this  report  will  also  be  modified  and 
submitted  for  publication. 
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The  Effect  of  Cementation  on  the  Elastic 
Properties  of  Granuleir  Material 

Jack  Dvorkin,  Gary  Mavko  and  Amos  Nur 


Abstract 

A  simple  analytical  model  is  developed  to  describe  the  mechanics  of  ce¬ 
mented  granular  material.  The  two-dimensional  problem  of  the  deforma¬ 
tion  of  an  arbitrarily-shaped  cement  layer  has  been  reduced  to  an  ordinary 
differential  equation  of  the  second  order.  A  simple  analytical  solution  is 
obtained  to  model  the  normal  and  shear  deformation  of  a  cement  layer  with 
straight  boundaries.  The  problem  of  the  normal  compression  of  two  con¬ 
tacting  circular  elastic  grains  and  an  elastic  cement  layer  between  them  has 
been  reduced  to  a  linear  integral  equation.  This  problem  has  been  solved 
both  for  two-dimensional  (cylindrical  grains)  and  three-dimensional  (spher¬ 
ical  grains)  cases.  The  elastic  modulus  of  such  a  combination  is  strongly 
increased  by  cementation  compared  to  the  case  of  a  classic  Hertzian  con¬ 
tact. 

Introduction 

The  acoustic  and  mechanical  characteristics  of  sedimentary  rocks  may  be 
strongly  affected  by  the  properties  and  structure  of  intergranular  bond  ma¬ 
terial.  The  deformational  pattern  of  cement  material  and  its  interaction 
with  elastic  grains  is  important  for  estimates  of  stress-strain  behavior  of 
granular  material,  as  well  as  for  its  failure  criteria.  Sedimentary  rocks  of 
interest  with  intergranular  cementation  span  a  wide  range  from  diagenetic 
sediments  and  sand-clay  mixtures  to  tar  sands. 

Numerous  publications  on  the  mechanics  of  granular  media  were  dis¬ 
cussed  by  Stoll  (1989).  The  theoretical  description  of  granular  material  has 
been  based  on  the  classical  solution  for  the  problem  of  a  normal  compression 
of  elastic  spheres  or  disks  by  Hertz  (Johnson,  1985)  as  well  aS  on  the  results 
regarding  the  oblique  compression  of  elastic  spheres  and  disks  (Mindlin, 


1949;  Walton,  1978).  All  these  models  consider  the  direct  contact  of  elas¬ 
tic  bodies.  The  dimension  of  a  contact  zone  changes  with  varying  external 
forces. 

The  situation  is  different  for  the  contact  of  initially  cemented  grains.  In 
this  case  the  dimension  of  a  contact  zone  between  two  bodies  is  predeter¬ 
mined  and  does  not  change  in  the  process  of  interaction.  Bruno  and  Nelson 
(1990)  examined  the  inelastic  mechanical  behavior  of  cemented  granular 
material  using  a  two-dimensional  discrete  element  procedure.  Intergranular 
interaction  through  cementation  was  modeled  using  a  linear  spring  scheme. 

In  this  paper  we  concentrate  on  the  description  of  cement  layer  defor¬ 
mation  as  well  as  on  the  interaction  of  elastic  grains  with  an  elastic  cement 
layer.  To  investigate  two-dimensional  cement  layer  deformation  when  the 
layer  is  thin  and  grains  are  much  stiifer  than  the  cement  we  have  derived  a 
simple  analytical  model.  The  approach  used  is  similar  to  those  commonly 
employed  in  approximate  theories  of  thin  plates  and  shells.  This  approach 
was  used  by  Matthewson  (1981)  to  model  an  axisymmetric  contact  on  thin 
compliant  coatings.  The  problem  of  a  two-dimensional  arbitrarily-shaped 
cement  layer  deformation  has  been  reduced  to  an  ordinary  differential  equa¬ 
tion  of  the  second  order.  A  simple  analytical  solution  has  been  derived  for 
the  case  of  a  cement  layer  with  straight  boundaries  under  compression  and 
shear. 

We  show  that  a  cement  layer  under  normal  compression  can  be  approx¬ 
imately  treated  as  an  Mastic  foundation”  between  two  grains.  Using  this 
result  we  reduce  the  two-dimensional  problem  of  interaction  between  two 
elastic  grains  and  an  elastic  cement  layer  to  a  linear  integral  equation.  Solv¬ 
ing  this  equation  we  are  able  to  estimate  the  elastic  modulus  of  cemented 
granular  material.  We  also  employ  the  Mastic  foundation”  approximation 
to  solve  the  problem  of  the  normal  compression  of  two  cemented  spherical 
elastic  grains.  The  presence  of  a  cement  layer  dramatically  increases  the 
elastic  modulus  compared  to  the  case  of  an  intergranular  Hertzian  contact. 

These  theoretical  models  can  be  used  to  predict  the  influence  of  cement 
content  and  location  on  the  properties  of  sediments. 

Deformation  of  a  Cement  Layer 
General  Equation 

To  examine  the  plane  deformation  of  an  elastic  cement  layer  we  introduce 
a  coordinate  system  {x,z)  in  the  cross-section  of  an  intergranular  boundary 


(Figure  1).  The  x  axis  is  directed  along  the  surface  of  the  lower  grain;  the  z 
axis  is  directed  upward.  The  problem  will  be  solved  in  a  linear  formulation, 
so  that  in  spite  of  the  possible  curvature  of  the  grain  surface  the  z  axis  can 
be  stiU  considered  as  being  directed  along  this  surface. 

Hook’s  law  for  the  plane  theory  of  elasticity  is: 


Oxx  =  (A  +  2a<)«x  + 

Oxx  —  A«i  +  (A  +  2^)wx, 

Oxt  =  /!(«»  +  tUr). 

The  equations  of  static  elasticity  are  the  following: 

dt^XT  .  ^^XZ  _  df^XZ  _  , 

17  ~  ’  17  17  ~ 


(1) 


(2) 


In  equations  (1)  and  (2)  <r„,  Ozz  and  a^x  are  normal  and  shear  stresses; 
u  and  w  are  displacements  in  the  z  and  z  directions;  lower  indexes  z  and  z 
denote  partial  derivatives;  A  and  fi  are  Lame’s  constants  of  the  cement. 

To  find  an  approximate  solution  for  the  problem  we  assume  that  the 
displacements  u  and  w  can  be  expressed  in  the  following  form: 


u(z,  z)  =  o(z)z  +  /3(z)z^  «;(z,  z)  =  ({x)z,  (3) 


where  the  functions  a,  0  and  (  are  to  be  found.  This  assumption  can  be 
iiuuitively  validated  by  the  deformation  pattern  of  initially  straight  fibers 
of  the  cement  (Fig.  2). 

We  also  assume  that  the  coordinate  system  (z,r)  is  connected  with  the 
surface  of  the  lower  grain,  thus  displacements  u(z,0)  =  u)(z,0)  =  0  at 
z  =  0.  The  surface  of  the  upper  grain  moves  relative  to  the  lower  grain,  and 
displacements  u  and  u'  are  given  along  this  surface: 


w  =  W(z);  u  =  U(x) 


at  z  =  /i(z),  where  h{x)  is  the  thickness  of  the  cement  layer. 

These  boundary  conditions  lead  to  the  following  relations  between  the 
functions  a,  0  and  c,  and  known  functions  kF(z),  U{x)  and  /i(z): 


f(z)  =  VF(z)//i(x); 
a{x)=V{x)lh{x)-0{x)h{x). 
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(4) 

(5) 


» 


I 


» 


I 


I 


Substituting  (3)  into  (1)  we  arrive  at  the  following  expressions  for  the 
stresses  in  the  cement  layer: 

=  (^  +  2/i)(a'2  +  +  At, 

Oxt  —  A(q^2  +  0*  z^)  +  (A  +  2/i)e,  (6) 

Oxx  =  M(a  +  2/32  +  iz). 

Integrating  the  first  equation  (2)  in  the  z  direction  from  0  to  h{x)  we 
arrive  at  the  following  relation: 

d 

Jo  +  = 

Substituting  expressions  (6)  into  this  relation  and  using  (5),  we  obtain 
the  following  equation  for  the  function  0(x): 

0"  +  Aix)0' +  Bix)0  =  C{x),  (7) 


where 


A  _  CLt/L.  D  _  ^fl.« 

yt  =  6* /k,  B  =  j|(, 

V  is  the  Poisson’s  ratio  of  the  cement. 

To  set  boundary  conditions  for  this  equation  we  will  introduce  normal 
force  N  in  the  x  direction: 

/V  =  r  a„dz  =  (A  +  2/i)(o'/iV2  +  /?'/iV3)  +  Ach. 

J  o 


This  force  is  assumed  to  be  zero  at  the  left  and  right  sides  of  the  cement 
layer;  A'  =  0  at  x  =  0  and  x  =  L,  where  L  is  the  length  of  the  layer  in  the 
X  direction. 

This  assumption  leads  to  the  following  boundary  conditions  for  the  equa¬ 
tion  (7): 


+  B— = -I- +  ■■■■  —  1 


(8) 


at  X  =  0  and  x  =  L. 

The  ordinacry  differential  equation  (7)  can  be  solved  numerically  using  the 
boundary  conditions  (8).  Once  the  functions  a  and  0  are  computed,  stresses 
Oxz,  Oxz  and  Oxx  can  be  found  from  relations  (6).  The  resulting  normal  and 
shear  forces  of  interaction  between  two  grains  are  to  be  computed  as  the 
integrals  of  stresses  in  the  x  direction  along  the  cement  layer. 
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Special  Case  -  Normal  Compression  of  a  Straight  Cement  Layer 

Considering  the  special  case  of  a  straight  cement  layer  under  compression 
we  have  the  following  conditions  for  the  functions  U{x),  H^(i)  and  h{x): 

f/  =  0;  VT  =  const-,  h  =  const. 


Equations  (7)  and  (8)  become: 


P"  -  ep  =  0; 


P'  =  C-,C  = 


del/ 

h^(l  -  v) 


at  I  =  0  and  x  =  L. 

Resolving  this  system  we  find 


_  C  —  e 

k  —  e“*^ 


If  the  length  of  the  cement  layer  significantly  exceeds  its  thickness  {L  >  h), 
we  have  the  following  estimates: 


kL  = 


L  6(1 -2u) 


>  1;  >  1;  <  1; 


1-1/ 


-kL 


Oxz  U=0«  -- 


^lW 


6i/2 


-fc(L-r)  _  g-fcxj 


h  ^(i-2u){i-uy 
Normal  stress  in  the  z  direction  Ozz  can  be  estimated  as  follows; 


»  (A  +  2fi)W/h. 

The  ratio  of  maximum  shear  stress  to  normal  stress  7  at  the  interface  of  the 
cement  layer  is  the  following: 

■  /3(l-2r^) 

'  i-i/y  2(1-1/)  ■ 

The  maximum  value  of  7  is  equal  0.47  at  1/  =  0.4.  Thus,  if  the  friction 
coefficient  between  the  cement  and  the  grain  surface  is  more  than  0.47,  the 
cement  layer  will  not  be  torn  from  the  grain  under  normal  loading. 
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Special  Case  -  Normal  Compression  and  Shear  of  a  Straight  Ce¬ 
ment  Layer 

In  this  case  we  consider  a  shear  displacement  superimposed  on  the  normal 
loading.  Namely,  we  put 

U  =  const  0;  W  =  const;  h  =  const. 

In  this  case  we  have  the  following  relations  between  a  and  (i:  '  . 

a  =  6-l3h,  a' =  -p'h. 

The  equation  and  boundary  condition  are  identical  to  the  ones  in  the 
previous  case  of  pure  normal  compression.  The  expression  for  the  normal 
stress  Ozi  does  not  change:  o,,  «  (A  +  2n)Wfh.  Shear  stress  Oxz  is:  = 

n{6  —  f3h)  at  2  =  0,  and  Oxz  =  +  /?h)  at  z  =  h.  In  this  case  the  maximum 
shear  stress  at  the  interface  “grain-cement”  is  approximately: 

U  W  I  ^ 

The  ratio  7  of  the  maximum  sheaur  stress  to  normal  stress  is: 

1-21/  U  I  ^  . 

■  2(1  -  V  (1  -  t'Kl  -  2uy' 

In  this  case  parameter  7  can  be  large  enough  depending  on  the  ratio  U/W. 
Thus,  the  separation  of  a  cement  layer  from  a  grain  surface  can  occur. 

Numerical  Examples  -  Tapered  Cement  Layer;  Layer  Between  Two 
Circles 

In  the  first  example  we  examined  the  deformation  of  a  tapered  cement  layer 
under  normal  compression  (Fig.  3).  The  length  of  the  layer  is  2  mm;  the 
minimal  thickness  of  the  layer  at  its  left  end  is  0.1  mm;  the  tangent  of 
its  tilt  is  0.2.  The  elastic  constants  of  the  cement  are:  Young’s  modulus 
E  =  2  ■  10^  Pa;  Poisson’s  ratio  1/  =  0.3.  The  normal  displacement  of  the 
cement  upper  surface  towards  its  lower  surface  is  10“^  mm.  Average  normal 
and  tangential  stress  distributions  along  the  cement  layer  are  given  in  Fig. 
3.  Actual  normal  stresses  (bold  curve)  are  compared  with  normal  stresses 
computed  using  an  approximate  formula  of  uniaxial  deformation:  Ozz  % 
(A  -I-  2n)W{x)lh{x)  (thin  curve).  These  two  curves  are  very  close. 


In  the  second  example  we  examined  the  deformation  of  a  cement  layer 
between  two  identical  circles  of  a  radius  1  mm  (Fig.  8,  a).  The  length  of 
the  layer  is  0.3  mm;  its  minimal  thickness  is  0.01  mm.  Elastic  constants  and 
compressional  displacement  are  the  same  as  in  the  first  example.  Average 
normal  and  tangential  stress  distributions  along  the  cement  layer  are  given 
in  Fig.  4.  In  this  case  actual  normal  stresses  (bold  curve)  are  also  compared 
with  normal  stresses  computed  using  the  uniaxial  deformation  approxima¬ 
tion  (thin  curve).  This  example  also  proves  the  accuracy  of  approximate 
uniaxial  deformation  formula. 

These  examples  show  that  noticeable  shear  stresses  may  develop  under 
normal  compression  in  a  cement  layer  of  changing  thickness.  Normal  stresses 
in  such  a  layer  can  be  approximately  calculated  using  the  assumption  that 
the  material  in  every  cross-section  of  the  layer  deforms  uniaxially:  <7,^  ss 
(^X  +  2(i)Wix)/h{x). 

Introducing  the  Elasticity  of  Grains  -  2D  Case 

The  solution  for  the  problem  of  two-dimensional  deformation  of  a  cement 
layer  between  two  grains  presented  above  can  be  incorporated  into  a  more 
general  problem  where  the  deformation  of  grains  is  taken  into  account.  The 
approximate  solution  of  such  a  problem  can  be  obtained  using  the  result  of 
the  previous  section: 


(T„  »  {X  +  2fi)W{x)/h{x).  (9) 

Assuming  as  well  that  shear  stresses  at  grain  surfaces  do  not  significantly 
influence  their  compressional  deformation  (Johnson,  1985)  we  conclude  that 
the  concept  of  “elastic  foundation”  {ibid.)  can  be  used  to  model  the  cement 
layer  compressional  deformation  between  two  deforming  grains  (Fig.  5). 

Examining  the  plane  deformation  of  an  elastic  body  on  a  cement  layer  - 
“elastic  foundation”  -  we  have  the  following  relation  between  the  displace¬ 
ment  of  the  cement  M^(i)  and  the  displacement  of  the  body  v(z)  (Fig.  6): 

t;(i)  =  €  +  »V(x),  (10) 

where  ^  is  the  rigid  translation  of  the  elastic  body.  Here  and  below  we 
consider  the  deformation  of  the  elastic  body  and  the  cement  layer  relative  to 
some  fixed  rigid  surface  (Fig.  6)  that  may  be  the  plane  of  symmetry  between 
two  identical  grains  (Fig.  8).  Assuming  that  stress  distribution  in  the  elastic 
body  is  close  to  one  in  a  half-plane  with  the  identical  surface  deformation 


(this  assumption  is  valid  provided  that  the  interface  “grain-cement”  is  much 
smaller  than  whole  grain  surface),  we  arrive  at  the  following  relation  between 
the  displacement  v(x)  and  normal  surface  stress  (Johnson,  1985): 

2(1  —  1/^)  /“ 

r(x)  = -  *  /  <r*j(s)  In  |x  —  a|ds -f  const,  (11) 

xIL/\  J—a 

where  is  Poisson's  ratio  and  Ei  is  Young’s  modulus  of  the  elastic  body; 
a  =  I./2  is  the  half-length  of  the  cement  layer;  the  origin  of  the  coordinate 
system  x  is  in  the  middle  of  the  cement  layer  (Fig.  7).  Equation  (11)  implies 
the  linearization  of  the  problem,  i.e.  the  assumption  that  the  contact  part 
of  the  grain  surface  can  be  approximated  by  the  interval  |x|  <  a  on  the  x 
axis. 

Substituting  (9)  and  (10)  into  (11)  we  arrive  at  the  following  integral 
equation: 

lF(x)  =  A  /  V  In  |x  —  s|ds  -f  const,  (12) 

J-a  n(s) 

2{l-u^,)(X  +  2n) 

tEi 

This  equation  can  be  solved  numerically  to  find  function  lY(x).  The  con¬ 
stant  in  the  right-hand  part  of  the  equation  will  be  calculated  using  the 
condition  of  a  given  integral  compressional  force  Fc  per  unit  length  of  cylin¬ 
drical  grains: 


Fc 


{X  +  2^i)W{x) 
h{x) 


Cementation  Between  Circular  Grains 

If  a  cement  layer  separates  two  identical  circular  grains  of  radius  R  (Fig.  8, 
a),  its  thickness  2h{x)  is  related  to  coordinate  x  as  follows: 

Hx)  SS  ho -i- 

where  2ho  is  the  minimal  separation  between  grains.  In  this  case  equation 
(12)  has  the  following  form: 

/“  VFfsl 

- -  In  |x  -  s|ds  -I-  const, 
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where  n  =  2holR.  A  very  important  special  case  is  the  one  where  grains 
have  a  point  contact  (Fig.  8,  b)  and,  thus,  ho  =  0.  Equation  (12)  in  this 
case  is: 

W{x)  =  2R\J‘ s\ds  + const.  (13) 

Here  we  cannot  neglect  the  deformation  of  grains  even  if  the  cement  is  very 
soft.  If  grains  are  absolutely  ri^d,  this  combination  will  not  deform  at  all. 
The  intriguing  question  is:  how  the  rigidity  of  this  combination  can  be  com¬ 
pared  to  the  rigidity  of  two-grain  combination  without  cement  between  them 
(Hertzian  contact).  The  principal  difference  between  the  example  under  con¬ 
sideration  and  the  classical  Hertzian  contact  is  that  in  the  first  case  grains 
will  keep  a  point  contact  under  increasing  loading,  whereas  the  Hertzian 
solution  implies  the  development  of  a  finite  contact  zone.  To  answer  this 
question  we  solved  equation  (13)  numerically  using  the  quadrature  method 
(Delves  and  Mohamed,  1985). 

We  examined  the  deformation  of  a  two-grain  combination  with  cemen¬ 
tation  between  them.  The  grains  of  radii  R  =  1  mm  have  a  point  contact. 
Young’s  modulus  of  grains  Ei  =2-10**;  Poisson’s  ratio  Ui  =  0.3. 

Stress  distributions  along  the  grain  surface  between  x  =  0  and  x  =  a 
(Fig.  8,  b)  are  given  in  Fig.  9.  In  this  example  (a  =  jR/3)  we  examined 
eight  different  cases  depending  on  the  ratio  of  cement  to  grain  stiffness  m  = 
E/Ei.  The  Poisson’s  ratio  of  the  cement  was  0.3.  The  shape  of  the  stress 
distribution  dramatically  changes  depending  on  the  value  of  m,  whereas  the 
integral  of  normal  stress  along  the  contact  interface  is  constant  and  equal  to 
the  compressive  force. 

The  shape  of  stress  distribution  between  cemented  grains  (Fig.  10,  bold 
curve)  is  completely  different  from  the  classical  Hertzian  case  (Fig.  10,  thin 
curve).  The  dimension  of  the  contact  interface  between  two  grains  is  also  dif¬ 
ferent:  in  the  Hertzian  case  this  dimension  increases  with  confining  pressure, 
whereas  the  contact  region  between  cemented  grains  remains  constant. 

To  calculate  the  deformation  and  the  stiffness  of  the  two-grain  combi¬ 
nation  we  used  the  solution  for  the  problem  of  deformation  of  an  elastic 
cylinder  subjected  to  external  stresses  (Novozhilov,  1961).  Stress  distribu¬ 
tion  at  the  grain  surface  was  used  to  find  the  deformation  of  the  system 
under  a  given  confining  pressure.  The  stiffnesses  of  cemented  and  nonce- 
mented  systems  are  compared  in  Fig.  11.  In  this  case  we  used  parameters: 
m  =  0.5;  a  =  R/30.  When  the  confining  pressure  is  small,  the  stiffness  of 
the  cemented  system  significantly  exceeds  the  stiffness  of  the  Hertzian  sys¬ 
tem.  The  stiffness  of  the  Hertzian  system  increases  with  confining  pressure; 
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80  does  the  contact  interface.  The  stiffness  of  the  Hertzian  system  reaches 
the  stiffness  of  the  cemented  system  at  the  confining  pressure  where  the  con¬ 
tact  interfaces  within  these  two  systems  become  approximately  equal.  The 
stiffness  of  the  cemented  system  does  not  depend  on  the  confining  pressure. 

The  stiffness  of  the  cemented  system  for  different  m  and  o  is  compared  to 
the  stiffness  of  the  Hertzian  system  (Fig.  12).  The  stiffness  of  the  cemented 
system  increases  with  the  length  of  the  cement  layer  and  with  the  relative 
stiffness  of  the  cement  m.  It  significantly  exceeds  the  stiffness  of  the  Hertzian 
system  in  the  range  of  confining  pressure  under  consideration. 


Introducing  the  Elasticity  of  Grains  -  3D  Case 

To  model  the  normal  compression  of  two  spherical  elastic  grains  interlaid 
with  an  elastic  cement  layer  we  employ  an  approach  similar  to  that  in  a  two- 
dimensional  case.  Examining  the  axisymmetrical  deformation  of  a  system  in 
a  meridional  plane  (Fig.  8)  we  assume  that  normal  stresses  in  a  cement  layer 
can  be  approximated  by  the  ‘‘elastic  foundation”  concept  and,  therefore, 
expressed  by  equation  (9):  as  (A  +  2fi)W{x)/h{x),  where  z  is  the  axis 

of  symmetry  and  the  x  axis  is  placed  in  the  meridional  plane.  We  also  use 
the  fact  that  shear  stresses  do  not  significantly  influence  the  compressional 
deformation  of  the  system  (Johnson,  1985).  In  the  case  under  consideration 
the  displacement  of  the  cement  TF(z)  and  normal  displacement  of  a  grain 
surface  v{x)  are  related  by  equation  (10):  i7(ar)  =  f  +  W(x),  where  (  is  the 
rigid  translation  of  a  grain. 

Following  the  classical  Hertzian  solution  we  assume  that  stress  distri¬ 
bution  in  the  elastic  grain  is  close  to  one  in  a  half-space  with  the  identical 
surface  deformation.  In  addition  we  linearize  the  problem  assuming  that  the 
contact  region  on  the  grain  surface  can  be  approximated  by  a  circle  |z|  <  a 
in  the  z  =  0  plane. 

These  assumptions  lead  to  the  following  relation  between  the  displace¬ 
ment  v{x)  and  normal  stress  at  the  grain  surface  (Timoshenko  and  Good- 
ier,  1970): 


t;(a:) 


2(\  —  i/?l  rxcoKfi+’y/a^-x^iin^  ifi 


(14) 


The  integration  in  this  equation  is  implemented  in  the  z  =  0  plane,  inside 
the  circle  |z|  <  a,  where  a  is  the  radius  of  the  cement  layer  (Fig.  13);  s  and 
ip  are  the  variables  of  integration. 
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The  problem  under  consideration  is  an  axisymmetricaJ  one,  thus  the 
function  Czz{s,<f>)  depends  only  on  the  r  coordinate  in  the  polar  coordinate 
system  (r,^)  in  the  z  =  0  plane.  This  argument  can  be  expressed  through 
arguments  s  and  ip  as  follows:  f  =  \/x^  +  —  2xa  cos  p.  This  means  that 

=  —  2x5 cos  p).  Substituting  this  expres¬ 

sion  into  equation  (14)  and  using  equations  (9)  and  (10)  we  arrive  at  the 
following  integral  equation  for  the  function  U^(z): 


fXCOZtfi- 


"f  ^o*— X*  •in*  ifi 


W(y/x^  +  s*  —  2zscos  (p)^^  ^ 

h{y/x^  +  —  2ZSCOS(,0)  ’ 


where  R  is  the  radius  of  the  grain;  A»  =  [4A(1  —  t'i)(A  +  2fi)]/{irEi).  Using 
the  fact  that  for  the  cement  layer  between  two  contacting  spherical  grains 
hix)  «  ^ ,  we  transform  this  equation  to  the  following  one: 


.  r  ,  /•*c<».^+v/a2-rWM>  Wiy/X^  +  S^  -2XS  COS  p)  ^ 

W{x)  =  -A,  dp  - 5-- — 5 — - ds  -  t. 

Jo  Jo  x^  +  —  2xs  cos  p 


(15) 

The  constant  ^  in  this  equation  can  be  calculated  using  the  condition  of  a 
given  integral  compressional  force  Fc  acting  on  the  grains: 


Fc  =  f  —<rzzix)2irxdx. 

Jo 


Equation  (15)  can  be  solved  numerically  similar  to  equation  (13).  In 
this  case  the  quadrature  method  is  applied  to  the  double  integral.  Once  the 
stress  distribution  along  the  grzun  surface  is  found,  the  deformation  and  the 
stiffness  of  the  two-grain  combination  can  be  calculated.  To  find  the  grain 
surface  displacements  we  used  the  solution  for  the  problem  of  deformation 
of  a  symmetrically  loaded  elastic  sphere  (Lur’e,  1964). 


Numerical  Examples  -  Cemented  Spherical  Grains 

We  examined  the  deformation  of  two  contacting  identical  cemented  elastic 
spheres  (Fig.  8.,  b).  All  constants  were  chosen  to  be  identical  to  those  in 
the  2D  case  above. 

Axisymmetrical  stress  distributions  at  the  grain  surface  differ  from  the 
Hertzian  stress  distribution;  they  are  qualitatively  similar  to  the  pattern 
observed  in  the  2D  case. 
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The  stiffness  of  cemented  system  is  compared  to  the  stiffness  of  a  Hertzian 
system  for  different  values  of  m  and  ajR  (Fig.  14).  The  stiffness  of  the  ce¬ 
mented  system  does  not  depend  on  the  confining  pressure.  It  exceeds  the 
stiffness  of  the  Hertzian  system  at  low  confining  pressure  even  when  the 
cement  is  soft  and  the  radius  of  the  cement  layer  is  small. 

Conclusions 

The  two-dimensional  problem  of  deformation  of  an  arbitrarily-shaped  ce¬ 
ment  layer  has  been  reduced  to  an  ordinary  differential  equation  of  the 
second  order.  A  simple  analytical  solution  has  been  derived  for  the  case 
of  a  cement  layer  with  straight  boundaries  under  compression  and  shear. 
A  cement  layer  under  normal  compression  can  be  treated  approximately  as 
an  “elastic  foundation”  between  two  grains.  Using  this  result  we  reduce  the 
problem  of  interaction  between  two  elastic  grains  and  an  elastic  cement  layer 
to  a  linear  integral  equation.  A  cement  layer  dramatically  increases  the  elas¬ 
tic  modulus  compared  to  the  case  of  intergranular  Hertzian  contact.  The 
shape  of  stress  distribution  between  cemented  grains  is  completely  different 
from  the  classical  Hertzian  case.  This  shape  dramatically  changes  depend¬ 
ing  on  the  ratio  of  cement  stiffness  to  grain  stiffness.  The  stiffness  of  the 
cemented  system  does  not  depend  on  the  confining  pressure.  The  stiffness 
of  the  cemented  system  increases  with  the  length  of  the  cement  layer  and 
with  the  relative  stiffness  of  the  cement.  It  significantly  exceeds  the  stiffness 
of  the  Hertzian  system  at  low  confining  pressure.  The  small  increase  of  the 
cementation  content  results  in  significant  growth  of  a  contact  zone  between 
two  contacting  grains  and,  thus,  dramatically  increases  the  stiffness  of  the 
system. 
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FIGURE  1:  Cement  Jayei  between  two  grains. 
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FIGURE  2;  Deformation  of  initially  straight  fibers  in  the  cement  tinder 
compression  and  shear. 
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FIGURE  3:  Average  norina]  and  shear  stress  distributions  in  the  tapered 
cement  layer  under  normal  compression. 


FIGURE  4;  Average  normal  and  shear  stress  distributions  in  the  cement 
layer  between  two  circles  under  normal  compression. 


FIGURE  7:  Th€  deformation  of  the  elastic  grain  assumed  to  be  dose  to 
the  deformation  of  an  elastic  half-plane. 
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FIGURE  8:  Cemcntalion  between  two  circular  grains  (cross-section  in  a 
two-dimensional  case;  meridional  plane  in  a  three-dimcnsiwiial  case). 
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FIGURE  10:  Norma]  stress  distributions  at  the  surface  of  two  contacting 
grains  with  the  cement  between  them  (bold  curve)  and  classical  Hertzian 
distribution  (thin  curve). 
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FIGURE  11:  The  stiffness  of  cemented 
pressure. 
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FIGURE  12:  The  stiffness  of  cemented  and  “Hertiian”  systems  vs.  con¬ 
fining  pressure  depending  on  the  ratio  of  cement  to  grain  stiffness  and  on 
the  length  of  the  cement  layer. 
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FIGURE  14:  The  stifioess  of  Gemented  spheres  systems  and  Hertzian  sys¬ 
tems  vs.  conhnJng  pressure  depending  on  the  ratio  of  cement  to  grain  stiff¬ 
ness  and  on  the  length  of  the  cement  layer. 
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Abstract 

We  present  a  theoretical  model  for  the  compaction  of  sediments  formed  by 
quartz  grains  separated  at  their  contacts  by  a  cementing  component.  Com¬ 
paction  occurs  due  to  lithostatic  pressure  that  results  in  increasing  contact 
stresses  with  burial.  As  the  stresses  at  the  contact  become  larger,  the  ce¬ 
menting  material  yields,  the  separation  between  the  grains  decreases,  and 
the  area  of  the  cemented  contact  increases.  This  results  in  an  increase  in 
P-  and  S-wave  velocity  with  depth.  The  rate  of  increase  depends  on  the 
strength  of  the  cement.  Similarly,  the  rate  of  porosity  reduction,  which  is 
dependent  on  the  reduction  of  separation  of  the  grains,  also  depends  on  the 
strength  of  the  cement.  Although  the  compaction  of  the  sediment  occurs 
due  to  plastic  yield  of  the  cementing  component,  the  material  as  a  whole 
behaves  elastically  in  response  to  seismic  wave  excitation.  Our  estimates  of 
the  seismic  properties  of  compacting  sediments  are  based  on  the  theory  of 
cemented  granular  material  (Dvorkin  et  al.,  1991)  that  is  used  to  compute 
the  contact  stiffness  of  cemented  grains.  Seismic  velocities  are  evaluated 
using  standard  formulas  that  relate  the  elastic  properties  of  granular  mate¬ 
rials  to  the  contact  stiffness  and  the  average  number  of  contacts  per  grain. 
The  results  indicate  that  the  velocities  steeply  increase  on  the  initial  stage 
of  compaction.  As  compaction  continues,  the  velocities  gradually  approach 
constant  values  that  correspond  to  the  direct  contact  of  the  grains.  If  the  ce¬ 
ment  is  extremely  soft  relative  to  the  grains,  a  velocity  maximum  occurs  just 
prior  to  the  pressure  at  which  the  grains  contact.  Our  theoretical  predictions 
agree  qualitatively  with  velocity-depth  profiles  within  oceanic  sediments. 

Introduction 

In  this  paper  we  develop  a  theory  for  the  compaction  of  granular  sediments 
due  to  increasing  lithostatic  pressure.  We  assume  that  rigid  grains  are  sep- 
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arated  by  a  cementation  material  that  is  being  pushed  out  of  grain  cont2u:ts 
as  confining  pressure  increases.  Using  simple  yield  criterion  we  find  the 
separation  between  grains  and  the  size  of  the  cemented  zone  at  the  grain 
contacts.  The  knowledge  of  the  geometry  of  the  grain  contacts  allows  us  to 
apply  our  theory  of  cemented  elastic  grains  (Dvorkin  et  al.,  1991).  Then  we 
use  formulas  that  relate  acoustic  velocities  in  random  sphere  packs  to  the 
normal  and  tangential  contact  stiffnesses  and  the  number  of  contacts  per 
grain  (Winkler,  1983). 

It  is  important  to  mention  that  in  spite  of  the  plastic  regime  of  the 
compaction  of  sediments,  we  assume  that  they  react  elastically  to  the  seismic 
wave  excitation. 

In  addition  we  give  approximate  formulas  for  velocities  in  cemented  sed¬ 
iments  obtained  under  the  assumption  that  grains  are  absolutely  rigid. 

Plastic  Compaction  of  Granular  Sediments 

We  assume  that  two  adjacent  spherical  gruns  are  separated  at  their  con¬ 
tact  by  a  given  constant  amount  of  cementing  material.  As  contact  pressure 
increases,  cement  is  being  pushed  out  of  the  contact  zone.  This  process  re¬ 
sults  in  smaller  intergranular  separation  and  increasing  area  of  the  cemented 
contact  zone  (Fig.  1). 

We  use  the  following  simple  speculation  in  order  to  find  the  geometry  of 
the  cemented  contact  zones  as  a  function  of  the  confining  pressure.  If  the 
confining  pressure  is  Pc  then  the  total  force  Ft  acting  upon  the  surface  of  a 
spherical  grain  of  radius  R  is  Ft  =  AttR^Pc-  This  force,  if  distributed  among 
C  contacts,  produces  a  forces  per  contact  Fc  =  Ft/C  =  4irR?Pc/C. 

The  average  normal  stress  <r  in  a  cement  disk  of  central  angle  0o  (Fig. 
2)  is; 

Fc  ^  4Pc 

""  itiRBoy  cel' 

This  formula  allows  us  to  relate  central  angle  Bo  to  the  confining  pressure 
and  the  yield  limit  of  the  cement  <r,: 
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Compaction  and  the  Contact  Geometry 

The  half  volume  of  the  cement  disk  of  central  angle  6  (Fig.  2)  is: 


n  =  +  I ). 

where  ho  is  the  half-distance  between  two  adjacent  grains. 

Assuming  that  the  volume  of  the  cement  at  the  contacts  does  not  change 
with  compaction  and  is  given  by  the  constamt  ratio  of  cement  to  grain  volume 
kv,  we  find  that  the  total  volume  of  cement  per  grain  is 


Vcem  —  • 


Thus  the  volume  Vc  can  be  calculated  as 


Vc 


= 


4itR^ 
3C  ■ 


Now  the  half-distance  between  two  adjacent  grains  ho  can  be  related  to 
the  radius  of  a  grain  R  as: 


*2-  m 

R  3c5jE  4  ■  ’  ' 

Equations  (1)  and  (2)  determine  the  geometry  of  the  cementation  at 
the  contact  depending  on  the  confining  pressure.  Given  this  geometry,  the 
normal  contact  stiffness  of  two  spherical  grains  can  be  calculated  using  our 
theory  of  intergranular  cementation  (Dvorkin  et  al.,  1991). 


Seismic  Velocities  and  Contact  Stiffnesses 


To  find  seismic  velocities  Vp  and  V,  as  functions  of  confining  pressure  we  use 
formulas  relating  Vp  and  V,  to  the  normal  and  tangential  contact  stiffnesses 
Sn  aJ>d  St  (Winkler,  1983): 


3C 

2Q-kRp 


(5n  +  |5,),  v;  = 


c 

207rJ?p 


(5.  +  ^S.). 


(3) 


where  p  is  the  density  of  the  grain  material.  These  formulas  were  derived 
for  the  random  packing  of  identical  elastic  spheres.  Stiffness  is  defined  as 


the  ratio  of  an  acting  force  to  the  displacement  it  produces.  A  convenient 
interpretation  of  (3)  in  terms  of  bulk  and  shear  moduli  K  and  G  is: 


Seismic  velocities  can  be  expressed  through  K  and  G  as: 


where  pr  is  the  density  of  the  rock. 

Given  the  Poisson's  ratio  of  the  rock  u,  the  shear  modulus  G  can  be 
expressed  through  the  bulk  modulus  K  as: 


G  = 


2(1  +  u)  • 


(6) 


Formulas  (4)  -  (6)  can  be  used  to  calculate  Vp  and  V,  given  the  normal 
stiffness  of  two  grains  and  the  Poisson’s  ratio. 


Seismic  Velocities  and  Confining  Pressure 

In  this  section  we  use  the  solution  for  the  problem  of  the  normal  compression 
of  cemented  circular  grains  (Dvorkin  et  al.,  1991)  in  order  to  find  Vp  and  V, 
versus  confining  pressure  Pc.  The  order  of  calculating  the  velocities  is  the 
following: 

equation  (1)  is  used  to  find  sis  a  function  of  Pci 

equation  (2)  is  used  to  find  the  half-distance  ho\ 

results  of  Dvorkin  et  al.  (1991)  are  used  to  find  normal  stiffness  Sni 

the  bulk  modulus  K  is  calculated  from  (4); 

a  certain  value  is  assigned  to  the  Poisson’s  ratio  v  of  the  rock; 

formulas  (6)  and  (5)  are  used  to  calculate  G,  and  Vp  and  V,. 

The  results  of  calculations  with  i/  =  0.3  are  presented  in  Fig.  3.  Different 
curves  Vp-Pc  and  V,-Pe  were  obtained  for  various  values  of  ratio  m  =  Ed  Eg, 
where  Ec  and  Eg  are  the  Young’s  moduli  of  the  grain  material  and  the 
cement.  We  also  varied  the  amount  of  cement  in  the  rock  (parameter  kv). 

Plots  of  Vp  and  V,  versus  Pc  for  the  case  where  m=0.07  ifc„=0.2  are 
presented  in  Fig.  3A.  Velocities  steeply  increase  on  the  initial  stage  of  com¬ 
paction.  As  compaction  continues,  the  velocities  gradually  approach  con¬ 
stant  values  that  correspond  to  the  direct  contact  of  grains.  Fig.  3B  gives 
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plots  of  Vp  and  V,  versus  Pc  for  m=0.07  Jbt,=0.1.  Velocitits  here  are  smaller 
than  in  the  first  examle,  that  is  due  to  reduced  amount  of  the  cementation. 

An  interesting  effect  of  peaks  in  P-  and  S- velocity  before  they  reach  their 
constant  values  is  shown  in  Fig.  3C.  In  this  case  the  cementation  is  very 
soft  compared  to  the  grain  material  (m=0.01}  and  ibv=0.1.  This  is  due  to 
the  fact  that  the  normal  stiffness  of  two  spherical  grains  separated  by  a  thin 
layer  of  a  soft  cement  can  be  higher  than  that  of  directly  contacting  grains. 

Our  theoretical  predictions  have  a  qualitative  agreement  with  velocity 
measurements  in  oceanic  sediments. 


Simplified  Formulas  for  5„  and  St 

Simple  formulas  for  evaluating  the  contact  stiffnesses  of  two  cemented  grains 
can  be  obtained  under  the  assumption  that  the  grains  are  absolutely  rigid. 
These  formulas  can  be  applied  only  if  the  cementation  is  very  soft  compared 
to  the  material  of  the  grains  and  if  grains  are  not  in  the  direct  contact. 

In  the  following  derivations  we  use  the  result  of  Dvorkin  et  al.  (1991) 
that  a  cement  layer  can  be  treated  approximately  as  an  elastic  foundation. 
Namely,  we  assume  that  a  normal  stress  in  the  cement  layer  (in  the  direction 
perpendicular  to  the  surface  of  the  grain)  is: 

<T  =  (A  +  2/i)W^/h,  (7) 

where  W  is  the  displacement  of  the  cement  normal  to  the  grain  surface;  h 
is  the  thickness  of  the  cement  layer  between  two  grains;  and  A  and  /i  are 
Lame’s  constants  of  the  cement. 

For  the  shear  stress  r  we  will  use  the  following  formula: 


r  =  ^,U/h,  (8) 

where  U  is  the  tangential  displacement  between  two  grains. 

The  thickness  of  the  cementation  between  two  circular  grains  is  (Fig.  2): 

/i(T)  =  2(ho  +  |i).  (9) 

From  (8)  and  (9)  we  find  that 


^  ^  hix)  *2  +  2Pho' 
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The  resulting  tangential  force  T  between  two  grains  produced  by  the 
displacement  U  is  for  the  2-D  case: 

T  =  2  J  T{x)dx  =  \j2R/ho  arctan(^^fi/ho). 

The  normal  force  produced  by  the  displacement  W  can  be  found  from 

(7): 


These  formulas  can  be  used  now  to  find  the  normal  and  tangential  stiff¬ 
nesses  between  two  grains  for  the  two-dimensional  case: 

St  =  T/(U/2)  =  2MV5ikarctan(tfo^/^), 

5„  =  Nf(W/2)  —  2(A  -i-  2/i)>/2r arctan(tfo^^/2), 
where  k  =  R/ho. 

For  the  3-D  case  of  two  cemented  spherical  grains  we  have: 

T=  f  2rxT(x)dx  =  irfiURla(l -i- 
Jo  2 


N  =  ir{X  +  2n)WRlnil  -H  ^Jb), 


St  =  2irfiUln{l  +  ^k\ 

(10) 

Sn  =  2w{X  +  2n)R]ii{l  +  ^k). 

(11) 

Velocities  Vp  and  V,  can  be  calculated  now  from  (3). 

Formulas  (10)  and  (11)  can  be  used  only  for  the  rough,  an  order  of 
magnitude  estimates  of  seismic  velocities. 

Conclusions 

We  developed  a  theoretical  model  for  the  compaction  of  sediments  formed 
by  ripd  grains  separated  at  their  contacts  by  a  cementing  component.  Com¬ 
paction  occurs  due  to  lithostatic  pressure  that  results  in  increasing  contact 
stresses  with  burial.  Our  estimates  of  the  seismic  properties  of  compacting 
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sediments  are  based  on  the  theory  of  cemented  granular  material  that  is 
used  to  compute  the  contact  stiffness  of  cemented  grains.  The  compaction 
of  sediments  occurs  as  a  plastic  process  but  they  react  elastically  to  the 
seismic  wave  excitation.  Velocities  Vp  and  V,  steeply  increase  on  the  ini¬ 
tial  stage  of  compaction.  As  compaction  continues,  the  velocities  gradually 
approach  constant  values  that  correspond  to  the  direct  contact  of  grains. 
An  interesting  effect  of  smooth  peaks  in  P-  and  S- velocity  before  they  reach 
their  constant  values  is  predicted  if  cementation  is  very  soft  compared  to  the 
grain  material  (e.g.  quartz  grains  in  sUt).  Our  theoretical  predictions  have 
a  qualitative  agreement  with  velocity  measurements  in  oceanic  sediments. 
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Figure  1  :  Deformation  of  the  cementation  between  two  spherical  grains 
due  to  increasing  confining  pressure. 
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Figure  2  :  Cementation  layer  on  a  spherical  grain. 


Figure  3  :  Vp  (bold  curves)  and  V,  (thin  curves)  versus  confining  pressure 
Pc.  A.  kv=0.2,  m=0.07;  B.  A;u=0.1,  m=0.07;  C.  fc„=0.1,  m=0.01. 
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t  Seismic  Velocities  in  Compacting  Sediments 

Jack  Dvorkin  and  Daniel  Moos  (both  at:  Department  of  Geophysics, 
Stanford  University,  Stanford,  CA  94305) 

We  present  a  theoretical  model  for  the  compaction  of  sediments 
I  formed  by  quartz  grains  separated  at  their  contacts  by  a  cementing 
component.  Compaction  occurs  due  to  lithostatic  pressure  that 
results  in  increasing  contact  stresses  with  burial.  As  the  stresses  at 
the  contact  become  larger,  the  cementing  material  yields,  the 
separation  between  the  grains  decreases,  and  the  area  of  the  cemented 
contact  increases.  This  results  in  an  increase  in  P-  and  S-wave 
>  velocity  with  depth.  The  rate  of  increase  depends  on  the  strength  of 
the  cement.  Similarly,  the  rate  of  porosity  reduction,  which  is 
dependent  on  the  reduction  of  separation  of  the  grains,  also  depends 
on  the  strength  of  the  cement. 

Although  ^e  compaction  of  the  sediment  occurs  due  to  plastic  yield 
of  the  cementing  component,  the  material  as  a  whole  behaves 
%  elastically  in  response  to  seismic  wave  excitation.  Our  estimates  of 
the  seismic  properties  of  such  a  material  are  based  on  the  theory  of  a 
cemented  granular  material  (Dvorkin,  et  al.,  1991)  that  is  used  to 
compute  the  contact  stiffness  of  cemented  grains.  Seismic  velocities 
are  evaluated  using  standard  formulas  that  relate  the  elastic  properties 
of  granular  materials  to  the  contact  stiffness  and  the  average  number 
I  of  contacts  per  grain. 

The  results  in^cate  that  the  velocities  increase  rapidly  on  the  initial 
stage  of  compaction.  As  compaction  continues,  the  velocities 
approach  constant  values  that  correspond  to  the  direct  contact  of  the 
grains.  If  the  cement  is  extremely  soft  relative  to  the  grains,  a 
velocity  maximum  occurs  just  prior  to  the  pressure  at  which  the 
^  grains  contact.  Our  theoretical  predictions  agree  qualitatively  with 
velocity-depth  profiles  within  oceanic  sediments.  However,  the 
importance  of  these  effects  depends  on  the  strengths  of  the  materials 
involved. 

Dvorkin,  J.,  Mavko,  G.  and  Nur,  A.,  1991,  The  Effect  of 
I  Cementation  on  the  Elastic  Properties  of  Granular  Material, 
Mechanics  of  Materials,  in  press. 
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Modeling  the  Granular  Structure  of 

Rocks 


Abstract 

The  interaction  of  rock  grains  with  friction  and  slippage  is  modeled  using 
an  analytical  solution  for  plain  elastic  deformation  of  cylinders.  The  model 
explains  the  difference  between  “static”  and  “dynamic”  elastic  moduli  of 
natural  rock.  The  numerical  experiments  performed  on  a  four-grain  unit 
predict  the  dependence  of  “static”  elastic  moduli  on  confining  pressure,  the 
friction  coefficient  between  grains  and  their  relative  initial  positions.  The 
“dynamic”  elastic  moduli  depend  mainly  on  the  relative  initial  positions  of 
grains  and  confining  pressure. 

Four-grain  units  with  varying  initial  configurations  and  friction  coeffi¬ 
cients  can  be  used  for  modeling  granular  material.  The  “static”  and  “dy¬ 
namic”  moduli  of  this  material  will  depend  on  the  statistical  distribution  of 
the  above-mentioned  properties. 

Introduction 

The  theoretical  description  of  granular  material  regarding  the  mechanical 
behavior  of  natural  rock  has  been  a  published  subject  for  more  than  30  years. 
In  many  publications  the  authors  employed  regular  packing  of  spherical 
grains  (Gassman,  1951,  Duffy  and  Mindlin,  1957,  Deresiewicz,  1958,  White, 
1965,  Walton,  1975). 

Digby  (1981)  investigated  the  behavior  of  identical,  randomly  stacked, 
spherical  particles.  Contacting  particles  are  initially  bonded  together  across 
small  areas.  The  solution  predicts  the  dependence  of  elastic  wave  speeds  on 
the  confining  pressure  and  adhesion  radii  of  contacting  particles. 

Walton  (1987)  derived  the  effective  incremental  elastic  moduL  of  random 
packing  of  identical  elastic  spheres.  The  spheres  are  assumed  to  be  either 
infinitely  rough  or  perfectly  smooth.  The  results  derived  are  applicable  to 
dense  packing. 


1 


Out  goal  is  to  examine  the  deformation  of  unconsolidated  rocks  under  dif¬ 
ferent  types  of  loading  conditions.  Using  the  exact  analytical  solution  for  the 
deformation  of  an  elastic  cylinder  under  arbitrary  external  forces  distributed 
along  its  perimeter  (Novozhilov,  1961)  we  model  the  interaction  of  circular 
grains  with  friction  and  slippage.  Employing  the  simple  combination  of  four 
identical  circular  grains  we  explain  the  observed  difference  between  “static” 
and  “dynamic”  elastic  moduli  of  unconsolidated  rocks  assuming  that  in  the 
later  case  no  slippage  occurs.  Considering  frictional  forces  we  explain  the 
hysteresis  effect  in  stress-strain  relationships  for  granular  materials. 


Deformation  of  Circular  Cylinder 

The  solution  of  the  plane  problem  for  a  circular  cylinder  (Novozhilov,  1961) 
is  based  on  the  general  theory  derived  for  plane  deformation  of  elastic  bodies. 
Here  we  use  the  following  relationship  between  displacements  u  and  v  and 
complex  functions  ^(r)  and  VC^)  iQ  &  complex  plane  s  =  x  -f  ty  (Fig.  1): 

«  +  it)  =  i^[A^(2)  -  z>pi(z)  -  ^>(2)],  (1) 

where  E  and  u  are  the  Young’s  modulus  and  Poisson’s  ratio  of  the  cylinder; 
A  =  Z  -  4u]  (,?/(r)  and  rl){2)  are  functions  conjugate  with  respect  to  ipf{2) 
and  t^’(r). 

The  complex  functions  v?(r)  and  ^(r)  can  be  found  using  boundary  con¬ 
ditions  for  normal  fp  and  tangential  fe  traction  on  the  cylinder  perimeter. 

The  functions  ip(2)  and  V’(-?)  are  regular  within  the  region  under  consid¬ 
eration,  and  can  be  represented  in  the  form  of  power  series 


s=(.-)  =  f  «»)  =  E  f •P'M  =  E 


k=1 


ksii 


k=l 


k-1 


(2) 


The  coefficients  o^c  and  Pk  are  connected  with  the  boundary  tractions  as 
follows: 


oi  =  -^1/2;  Q*  =  Ak,  k>  \\  ^k  —  A-k  -  (^  +  2)Ak+2- 


^  ^  (/p  +  ifs)  €xp[-i(l:  -  1  )^]  de. 

Here  R  and  6  are  the  modulus  and  argument  in  the  complex  z  plane. 
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Conditions  at  Grain  Contact 


To  describe  a  aormal  stress  distribution  along  contact  surfaces  we  use  the 
Hertzian  approach  (Johnson,  1984).  The  contact  stripe  between  two  cylin¬ 
ders  pressed  against  each  other  by  a  force  N  has  a  width  2a,  where 

- ’  (3) 

and  Rc  is  the  cylinder  radius. 

Pressure  distribution  along  the  contact  surface  can  be  presented  in  the 
following  form: 

/.({) = (4) 

w’here  ^  is  the  local  coordinate  system  (see  Fig.  2). 

To  model  the  distribution  of  a  tangentional  traction  along  the  contact 
surface  of  two  cylinders  when  slippage  occurs  we  use  the  Coulomb  law  of 
friction: 


/«  =  «/<>» 

where  k  is  the  coefficient  of  sh'ding  friction. 

The  modeling  of  tangential  traction  distribution  in  the  absence  of  sliding 
is  more  complicated.  Mindlin's  (1954)  solution  of  this  problem  predicts 
that  if  two  spheres  are  pressed  together  by  a  constant  normal  force  and 
then  tangential  forces  are  appLed,  the  pattern  of  normal  traction  remains 
unchanged  while  the  resulting  distribution  of  tangential  tractions  tends  to 
infinity  at  the  periphery  of  the  contact  area.  Since  infinite  stresses  are 
not  physically  possible,  stress  relief  in  the  region  of  stress  concentration  is 
allowed.  This  effect  gives  rise  to  a  tangential  traction  distribution  that  is 
not  proportional  to  the  normal  traction. 

Walton  (1987)  obtained  part  of  the  solution  corresponding  to  shear  using 
the  assumption  of  perfect  adhesion:  that  is,  no  relative  slip  between  the 
two  spheres  is  allowed  over  the  contact  area,  while  the  final  state  is  being 
attained.  The  resultant  configuration  and  the  occurence  of  slippage  are 
strongly  dependent  on  whether  the  two  spheres  are  first  compressed  and 
then  sheared  or  whether  the  two  motions  are  occuring  simultaneously.  In 
this  case  the  tangential  traction  distribution  has  a  form  similar  to  the  one 
given  by  £qn.  4.  In  this  case  if  slippage  occurs,  it  will  be  in  the  form  of 
sliding  over  the  whole  of  the  contact  area. 
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In  our  model  we  use  Walton’s  solution  for  the  shear  traction  distribution: 


/,(i)  =  (5) 

where  T  is  the  resulting  tangential  force  applied  to  the  contact  area. 

The  dimensions  of  the  contact  area  between  two  cylinders  is  usually 
very  small  compared  to  cylinder  size.  This  allows  us  to  use  in  this  case 
Saint- Venant’s  principle,  according  to  which  statically  equivalent  systems 
of  forces  distributed  over  a  small  part  of  the  boundary  of  an  elastic  body 
produce  in  the  body  systems  of  stresses  which  differ  from  each  other  only 
in  the  immediate  vicinity  of  the  area  of  loading.  However,  at  sufficient 
distances  from  this  area,  the  stresses  aroused  by  one  or  the  other  system  of 
forces  are  practically  identical.  Thus,  under  certain  conditions  the  form  of 
shear  traction  distribution  is  not  really  important  to  compute  the  resulting 
deformation  of  cylindrical  grains. 

Four  Cylinders  System  -  Active  and  Reactive  Forces 

To  understand  the  peculiarities  of  grain  interaction  we  examine  a  simple 
combination  of  four  identical  cylinders  (Fig.  3A)  subjected  to  the  internal 
loadings  Pi  and  P^. 

The  equilibrium  condition  of  four  cylinders  system  can  be  presented  in 
the  following  form: 

Pi  =  2{N cosq  +  T sin o),  P2  =  2(A’ sin o  —  Tcos q),  (6) 

where  N  and  T  are  the  normal  and  tangential  reactive  forces  acting  between 
two  cylinders  (Fig.  3B);  o  is  an  angle  defining  the  relative  position  of  the 
cylinders. 

When  slippage  occurs  the  relationship  between  N  and  T  is  defined  by 
the  Coulomb  law  of  friction:  |  T  |=  kN.  When  there  is  no  relative  motion 
between  two  cylinders  the  reactive  forces  are  connected  by  the  following 
formula:  |  T  |<  k.N.  Making  use  of  Eqn.  7  we  arrive  now  at  the  following 
relation  between  the  internal  forces  Pi  and  P2  provided  there  is  no  slippage: 

(tana  -  k)/(H- « tan q).  (7) 

Once  the  slip  occurs  this  relation  turns  into  the  following  strict  equation: 

Fa  =  Fi(tana  —  ic)/(H- «tano).  (8) 
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Considering  P2  in  E^n.  9  as  the  function  of  tana  and  taking  the  first  deriva¬ 
tive  with  respect  to  this  argument  we  can  arrive  at  a  positive  expression. 
Thus,  to  achieve  an  equilibrium  of  the  four  cylinders  the  side  force  P2  must 
increase  with  increasing  a.  It  means  that  if  P2  is  constant,  and  P\  exceeds 
a  limiting  value  dictated  by  Eqn.  9  the  system  cannot  stay  in  balance.  In 
this  case  the  system  collapses  and  the  angle  a  changes  from  its  initial  value 
to  its  upper  limit  of  60°. 

The  system  collapse  can  be  prevented  provided  that  displacement  bound¬ 
ary  conditions  are  used.  We  can  simulate  the  axial  loading  of  the  system 
increasing  the  vertical  force  Pi  and  assuming  that  the  system  is  limited  in 
its  expansion,  so  that  the  points  S\  and  S2  of  the  side  cylinders  (Fig.  4) 
do  not  suffer  any  displacements  that  can  increase  the  system  width.  Thus, 
we  arrive  at  the  uniaxial  deformation  of  the  system.  This  computational 
experiment  can  be  used  to  find  the  elastic  modulus  M  =  A  -f  2/x,  where  A 
and  fi  denote  the  Lame  elastic  constants  of  the  cylinders  material,  and,  thus, 
to  estimate  P-wave  velocity. 

To  perform  this  experiment  we  applied  at  the  be^nning  a  certain  con¬ 
fining  pressure,  and  then  gradually  increased  the  vertical  force  Pi.  Once 
slippage  occured  Eqn.  9  was  used  to  calculate  the  force  Pj  during  the 
loading  process.  The  angle  o  was  changed  to  meet  the  conditions  of  zero 
displacements  for  the  points  5i  and  5}. 

Unloading  the  system  we  gradually  decreased  the  vertical  force  Pi  check¬ 
ing  the  slippage  condition  (Eqn.  6).  To  find  the  force  P2  during  the  slippage, 
in  the  case  when  the  tangential  force  T  is  directed  to  prevent  the  cylinder's 
upward  motion,  we  used  the  following  equation  instead  of  Eqn.  9: 

Pj  =  Pi(tana -1- k)/(1  —  Ktano).  (9) 

Eqns:  1  -  9  can  be  used  to  compute  the  displacements  of  the  system 
during  loading  and  unloading. 

Numerical  Examples 

The  computations  were  performed  for  the  system  of  glass  cylinders  (E=62.27 
GPa,  i/=0.2466);  the  cylinder  radius  was  0.0375  mm.  These  parameters  were 
used  by  Marion  (1990)  in  his  experiments  conducted  on  ^ass  beads.  Friction 
coefficient,  initial  confining  pressure  (P^)  and  initial  geometry  of  the  system 
were  changed  in  our  computational  experiments  to  find  the  influence  of  these 
parameters  on  the  elastic  constants. 
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Four  Cylinders  System  -  Uniaxial  Deformation 

The  stress-strain  relationship  for  the  four  cylinders  systena  during  loading 
and  unloading  is  presented  in  Fig.  5.  In  this  case  the  initial  parameters 
values  were  the  following:  a  =  57®,  #c  =  0.3,  Pc  =  0.  The  initial  value  of 
angle  a  was  about  its  upper  limit  of  60®.  A  vertical  stress  a  was  calculated 
as  a  ratio  of  the  vertical  force  Px  to  the  width  of  the  system.  A  displacement 
t  was  found  as  the  ratio  of  a  vertical  displacement  and  the  initial  vertical 
size  of  the  system. 

The  stress-strain  curve  presented  in  Fig.  5  is  similar  to  the  theoretical 
curves  obtained  for  the  uniaxial  strain  loading  of  spheres  in  a  face-centered 
cubic  array,  and  the  results  of  experiments  performed  on  Ottawa  sand  (Stoll, 
1989).  The  hysteresis  effect  appears  due  to  the  change  of  friction  force 
direction  at  the  contact  area  between  adjacent  cylinders.  The  initial  pattern 
of  the  system  deformation  under  loading  condition  is  elastic.  After  sliding 
begins  the  system  deforms  plastically  with  the  angle  a  increasing.  In  the 
beginning  stage  of  unloading  the  upper  cylinder  is  stuck  between  the  side 
cylinders  and  the  system  deforms  elastically.  Ultimately,  when  sliding  in  the 
upward  direction  begins,  the  system  deforms  plastically. 

The  relationship  between  modulus  M  and  vertical  stress  a  is  presented 
in  Fig.  6  for  the  case  a  =  57®,  k  =  0.3  for  two  different  values  of  the  initial 
confining  pressure:  0  and  10  MPa.  The  curve  A  corresponds  to  the  curve  in 
Fig.  5.  Vertical  lines  connecting  the  bottom  and  top  branches  of  the  curve 
A  represent  the  change  of  M  when  unloading  is  initiated  at  certain  points 
of  loading  path  and  the  system  deforms  elasticaUy.  These  increased  values 
of  M  are  very  close  to  those  computed  during  the  monotonous  unloading  of 
the  four  cylinders  system  (the  top  branch  of  the  curve  A).  The  difference 
between  the  bottom  and  the  top  branches  of  the  curve  A  corresponds  to  a 
difference  between  “static”  and  “dynamic”  elastic  moduli  of  a  grain  material. 
In  this  example  the  ratio  of  the  “dynamic”  M  to  the  “static”  one  has  the 
order  of  5. 

The  curve  B  in  Fig.  6  corresponds  to  the  case  when  the  initial  confining 
pressure  is  10  MPa.  The  system  becomes  less  compliant  than  in  the  case 
when  Pc  =  0.  The  last  part  of  the  unloading  path  (the  top  branch  of  the 
curve  B)  is  close  to  an  elastic  pattern  unlike  the  last  part  of  the  curve  A 
when  intensive  plastic  deformation  occurs. 

The  computational  uniaxial  compression  experiments  were  conducted  on 
the  four  cylinders  system  for  different  angles  a,  friction  coefficients  k  and 
initiad  confining  pressure  Pg. 
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The  loading-nnloading  curves  for  o  =  57®,  Pc  =  0  and  k  =  0.6  and  1.2 
are  given  in  Fig.  7. 

The  increase  of  the  friction  coefficient  k  results  in  increased  "static” 
modulus  (curve  B,  the  bottom  branch)  and  decreased  "dynamic”  modulus 
(curve  B,  the  top  bratnch).  The  latter  effect  is  not  as  strong  as  the  former 
one.  The  origin  of  observed  difference  between  two  unloading  paths  (curves 
A  and  B,  the  top  branches)  is  the  loading  history. 

The  case  of  a  =  45®  and  ^  presented  in  Fig.  8.  The  values 
of  friction  coefficient  were  0.6,  0.9,  and  1.2.  In  the  last  case  the  system 
deformed  elastically  without  a  difference  between  loading  and  unloading 
patterns  (curve  C).  Loading  curves  correspond  to  plastic  deformation  in  the 
cases  when  k  =  0.6,  and  k  =  0.9.  The  "static”  modulus  M  increases  with 
increasing  k. 

For  a  s=  30®  (the  lower  limit  of  this  parameter)  the  system  deforms 
elastically  when  the  friction  coefficient  k  is  as  small  as  0.6  (curve  B,  Fig. 
9).  The  system  deforms  elastically  even  at  k  s  0.3  provided  that  initial 
confining  pressure  is  applied  (curve  C,  Fig.  9). 

If  K  =  0.3  and  P^  =  0  the  system  shows  plastic  behavior  when  loaded 
(the  bottom  branch  of  curve  A,  Fig.  9).  The  “dynamic”  moduli  calculated 
along  the  loading  path  are  very  close  to  those  obtained  from  the  elastic 
unloading  path  (the  top  branch  of  curve  A,  Fig.  9). 

The  above  examples  indicate  the  increasing  stiffness  of  the  system  when 
angle  a  decreases.  The  difference  between  the  “static”  modulus  (Jl/«)  and 
the  “dynamic”  one  (Md)  also  decreases  ,  so  that  the  ratio  Md/M,  decreases 
from  the  order  of  5  for  o  =  57®  to  the  order  of  1.5  for  o  =  30®. 

The  change  of  angle  a  along  the  plastic  deformation  curves  is  very  small 
and  has  the  order  of  0.5°. 

Four  Cylinders  System  >  Bulk  Modulus 

To  calculate  a  bulk  modulus  K  of  the  four  cylinders  system  we  apply  equal 
vertical  and  horizontal  forces:  Pi  =  P2.  Making  use  of  this  condition  and 
Eqns.  6  we  arrive  at  the  following  expression  for  the  ratio  of  tangential  to 
normal  force  at  the  contact  zone: 

7^  __  tano  —  1 
iV  ~  tan  Q  +  1  ’ 

This  relation  gives  the  following  range  of  T/N  variation  when  an^e  q 
changes  between  30°  and  60°;  —0.268  <  T/A  <  0.268.  The  realistic  values 
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of  inction  coeflicieDt  k  do  not  fall  into  this  interval.  Thus,  the  deformation 
of  the  system  when  equal  vertical  and  horizontal  forces  are  applied,  follows 
an  elastic  pattern. 

The  relationships  between  the  bulk  modulus  K  and  confining  pressure  Pc 
for  different  a  are  presented  in  Fig.  10.  The  values  of  K  increase  significantly 
when  a  decreases. 

Reducing  System^s  Stiffness 

To  reduce  the  stifihess  of  four  cylinders  system  we  assume  that  it  is  sur¬ 
rounded  by  rows  of  identical  cylinders  with  zero  displacements  of  the  right 
and  left  external  points  (Fig.  11). 

The  results  of  uniaxial  strain  experiments  for  such  a  system  are  presented 
in  Fig.  12.  In  this  case  the  four  cylinders  system  was  surrounded  by  twelve 
cylinders  on  both  sides.  The  values  of  angle  o  and  friction  coefficient  were 
the  following:  a  =  45®,  k  =  0.6. 

The  reduced  stiffiiess  of  the  system  strongly  affects  the  “static”  modulus 
M  (the  loading  curve.  Fig.  12).  The  reduction  of  this  parameter  has  the 
order  of  5  for  the  example  under  consideration. 

Discussion 

Making  use  of  an  exact  analytical  solution  for  the  plane  deformation  of  an 
elastic  cylinder  and  conducting  computational  experiments  on  the  simple 
combinations  of  cylinders  we  observed  the  difference  between  “static”  and 
“dynamic”  elastic  moduli  for  grain  material.  This  approach  allows  the  em¬ 
ployment  of  more  complicated  and  realistic  combinations  of  cylinders  for 
modeling  granular  material. 

The  simplest  approach  to  model  a  granular  aggregate  is  to  combine  four- 
grain  units  with  different  friction  coefficients  and  initial  configurations  in 
rows  and  columns.  The  distribution  of  mentioned  parameters  will  depend 
on  actual  properties  of  unconsolidated  rocks.  The  numerical  experiments 
conducted  on  such  models  result  in  relationships  similar  to  those  presented 
in  Figs.  5-10.  The  values  of  elastic  moduli  for  the  aggregates  are  equal  to 
weighted  averages  of  these  parameters  for  single  four-grain  units. 

Conclusions 

Considering  the  system  of  four  elastic  cylinders  interacting  with  friction  we 
discover  the  difference  between  “static”  amd  “dynamic”  elastic  moduU  when 
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slippage  occurs. 

The  “static”  modulus  of  the  system  depends  stron^y  on  initial  configu¬ 
ration,  friction  coefficient  and  confining  pressure.  The  “dynamic”  modulus 
depends  mainly  on  initial  configuration  and  confining  pressure. 

The  system’s  stiffiiess  can  be  significantly  reduced  due  to  mild  bo\mdary 
conditions  (four  cylinders  surrounded  by  rows  of  cylinders). 

The  elastic  properties  of  aggregates  of  four-grain  units  having  different 
initial  configurations  and  friction  coefficients  are  similar  to  those  of  single 
units.  The  values  of  elastic  moduli  are  equal  to  weighted  averages  of  these 
moduli  for  four-grain  units. 
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Figure  3  :  A.  Four  cylinders  combination.  B.  Normal  and  tangential 
tive  forces  between  two  cylinders. 


Figure  4  :  Four  cylinders  system  -  imiaxial  deformation. 


Figure  5  :  Stress-strain  relationship  for  the  uniaxial  deformation  of  four 
cylinder  system;  a  =  1,  k  =  0.3,  Pc  =  0- 


Figure  7  :  M  —  P  relationsiiip  for  the  four  cylinders  system’  a  =  57®  P  = 
♦  0,  K  =  0.6  and  1.2.  ’  ’  ‘ 


Figure  8  :  M  -  P  relationship  for  the  four  cylinders  system;  a  =  45°, 
Pj  =s  0,  K  =  0.6,  0.9,  and  1.2. 


1)  :  n  =  0,  K  =  O.C 
C:  re=  lOMPo,  «  =  0.3 


10  20  P,MPa 


Figure  B  :  M  -  P  relationship  for  the  four  cylinders  system;  o  =  3( 
Pc  =  0,  and  10  MPa,  k  =  0.3,  and  0.6. 
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Figure  10  ;  K  versus  for  a=30“,  45®,  and  60®. 


Figure  1 1  :  The  system  of  reduced  stifcess. 


Figure  12  :  M  -  P  relationship  for  the  system  of  reduced  stiffness;  four 
cylinders  system  is  surrounded  by  12  cylinders  on  both  sides;  q  :=  45®, 

K  =  0.6. 


The  Mechanics  of  Hollow  Grains 


Abstract 

We  develop  a  theory  of  the  two-dimensional  deformation  of  hollow  grains 
under  external  loading,  modeling  them  as  cylindrical  shells  of  a  closed  cross- 
section.  This  theory  allows  us  to  describe  the  deformation  of  arbitrarily- 
shaped  grains  with  walls  of  varying  thickness,  and  external  and  internal 
fractures.  Fractures  dramatically  reduce  the  stiffness  of  the  grains  and  may 
result  in  the  collapse  and  compaction  of  hollow  grains  under  lithostatic  pres¬ 
sure.  Calculating  the  normal  and  shear  stiffness  of  hollow  grains,  one  will 
be  able  to  determine  effective  elastic  properties  of  their  aggregate  using  the 
expressions  developed  for  a  random  packing  of  spheres.  The  theory  pre¬ 
sented  in  this  paper  can  be  directly  applied  to  estimating  the  acoustical  and 
mechanical  properties  of  basaltic  pillows  and  pelagic  sediments. 

Introduction 

The  acoustic  and  mechanical  characteristics  of  oceanic  crust  and  seabottom 
sediments  may  be  strongly  affected  by  the  properties  of  individual  grains 
that  compose  the  material.  Almost  all  theoretical  models  for  the  deforma¬ 
tion  of  a  granular  material  treat  an  individual  grain  as  an  elastic  element 
without  voids,  typically  as  a  disk  (2D)  or  a  sphere  (3D).  Recently,  Bruno 
and  Nelson  (1990)  examined  the  effect  of  intergranular  microfractures  on  the 
inelastic  behavior  of  sedimentary  rock.  The  deformation  of  the  aggregate 
of  arbitrarily-shaped  grains  can,  in  principle,  be  modeled  using  the  numer¬ 
ical  simulation  of  the  deformation  of  an  individual  grain  combined  with  a 
numerical  model  for  granular  assemblies  (Cundall  and  Strack,  1979). 

However,  there  are  at  least  two  cases  where  grains  with  interior  voids 
have  to  be  considered:  (1)  basaltic  pillows  that  are  present  in  the  oceanic 
crust,  and  (2)  pelagic  sediments  formed  by  hollow  sheUs.  In  this  paper  we 
concentrate  on  the  theoretical  description  of  the  deformation  of  a  hollow 
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arbitrarily-shaped  grain  under  given  loading.  We  model  a  grain  as  a  cylin¬ 
drical  shell  of  arbitrarily-shaped  middle  surface  and  varying  thickness.  This 
approach  allows  us  to  introduce  external  and  internal  cracks  in  the  walls  of 
the  grain  as  well  as  to  examine  the  effect  of  differently-shaped  internal  holes. 
The  problem  is  reduced  to  a  system  of  ordinary  differential  equations  that 
requires  numerical  solution.  We  show  how  the  shape  and  fracturing  affect 
the  stiffness  of  grains. 

This  is  a  two-dimensional  theoretical  approach,  but  it  can  quantitatively 
predict  the  characteristics  of  elongated  basaltic  pillows,  and  can  give  quali¬ 
tative  estimates  of  behavior  of  topologically  complicated  pelagic  sediments. 

Deformation  of  Cylindrical  Shells 

In  this  section  we  derive  a  general  theory  for  the  2-D  deformation  of  cylindri¬ 
cal  shells  using  assumptions  common  for  the  theory  of  thin  shells.  Typically, 
a  shell  is  considered  to  be  thin  if  the  ratio  of  its  thickness  2h  to  the  radius  of 
curvature  of  its  middle  surface  R  is  less  than  .05.  This  condition  guarantees 
the  accuracy  of  results  within  a  5%  range.  However,  the  theory  of  thin  shells 
can  be  also  applied  (with  decreased  accuracy)  for  thick  shells  with  the  ratio 
2h/R  having  the  magnitude  up  to  1/3  (Grigoluk  and  Kabanov,  1978).  Thus 
the  theory  of  deformation  of  hollow  gfains  based  on  the  theory  of  shells  is 
accurate  enough  to  quantitatively  estimate  the  acoustical  and  mechanical 
properties  of  basaltic  pillows  and  pelagic  sediments. 

Displacements  and  Forces  in  Cylindrical  Shells 

To  examine  the  2-D  deformation  of  an  infinitely  long  cylindrical  shell  (Fig. 
lA),  we  introduce  a  local  cylindrical  coordinate  system  (ff,r)  in  the  normal 
cross-section  of  the  shell  (Fig.  IB).  The  center  of  the  (ff ,  r)  coordinate  system 
is  placed  in  the  local  center  of  curvature  of  the  middle  surface  of  the  shell. 
We  also  introduce  a  local  rectangular  coordinate  system  (i,  z)  with  the  origin 
in  the  middle  surface  of  the  shell,  the  x  axis  tangent  to  this  surface,  and  the 
z  axis  directed  along  the  radius  of  curvature.  These  two  coordinate  systems 
are  related  to  each  other  as: 


X  =  R0;  z  =  r~  R.  (1) 

The  displacements  of  the  shell  are  u(x,z)  in  the  x  direction  and  tn(i,a) 
in  the  z  direction.  We  assume  the  following  form  of  functional  dependence 
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of  these  displacements  on  the  x  and  y  arguments: 


=  u®(x)  +  z4>{x)\  ib{x,z)  =  (2) 

This  assumption  implies  that  a  straight  fiber  initially  perpendicular  to  the 
middle  surface  keeps  its  straight  shape  after  deformation  and  tilts  at  an^e 
^  to  its  initial  direction  (Fig.  2A),  and  the  thickness  of  the  shell  does  not 
change. 

We  examine  the  case  of  a  plane  deformation  of  the  sheU:  tyy  =  0,  where 
Cyy  is  a  normal  deformation  in  the  y  direction  perpendicular  to  the  plane  of 
deformation.  This  condition  and  Hook’s  law  yield: 


_  ^.grr-t-ggg  _  Q 

E  E 


(3) 


where  E  is  Young’s  modulus;  v  is  Poisson’s  ratio;  Oyy,  Orr  and  obb  are 
components  of  a  stress  tensor  in  the  shell. 

Commonly,  in  the  theory  of  plates  and  shells  it  is  assumed  that  the  nor¬ 
mal  stress  in  the  direction  perpendicular  to  the  middle  surface  is  small  com¬ 
pared  to  the  other  normal  stresses  (Timoshenko  and  Woinowsky-Krieger, 
1959):  Ott  +  obb  «  obb-  This  assumption  and  (3)  yield: 


Oyy  =  vobb.  (4) 

Using  (4)  and  Hook’s  law  we  arrive  at  the  following  relations  between 
the  components  of  the  stress  and  deformation  tensors: 


The  components  of  the  deformation  tensor  can  be  related  to  the  dis¬ 
placements  ii  and  w,  and  functions  u*’,  vf  and  in  the  (5,r)  cylindrical 
coordinate  system  as: 

l.au'’  ,  d<t>  , 

««  = -(gj +  «.)=;(  ■p-  +  »5j  +  «>), 


1 dw  .  du 


^  ,  %JU  V  X  UW  V  QfZ 

“  719  ^l9~7'^  7W  +  ^ 

Using  (3)  and  assuming  that  r  »  R  in  the  shell,  we  transform  these 
equations  to: 

du°  d4>  vf  _  dw**  d>z 
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We  introduce  integral  normal  N  and  shear  Q  forces  in  the  shell,  and 
bending  moment  M  (Fig.  2B)  created  by  normal  stresses  as: 

N  =  f  tTgffdz;  Q  =  f  (Trfidz;  M  =  f  aggzdz.  (7) 

J-h  J-h  J-h 

Using  these  definitions  and  relations  (5)  and  (6),  we  express  the  integral 
forces  and  the  bending  moment  through  functions  v'*,  and  <l>: 


^  Eh  ,dw^  ^  ii",  „  2Eh  ,du^  tc»,  2Eh^  dd> 

«  =  TTr<-5r+*-:R  >■  w 


It  is  important  to  mention  that  the  x  coordinate  can  be  continuously 
counted  along  the  middle  surface  of  the  shell  by  varying  x  between  0  and  X, 
where  X  is  the  shell’s  perimeter.  This  approach  will  not  change  the  above- 
derived  formulas  as  they  include  only  the  differential  dx.  Parameters  R  and 
h  have  to  be  considered  as  functions  of  x. 


Equations  of  ShelPs  Deformation 

The  equations  for  the  deformation  of  the  shell  under  external  forces  will  be 
derived  from  the  equations  of  static  dasticity: 


d{T(j 

rr  )  . 


=  0; 


(9) 


d(rcrg)  ,  d<Tge  .  _  „ 

Integrating  (9)  from  —h  to  h  in  the  z  direction  and  using  (1)  and  (7)  we 
have: 

R(<J..U  -  <frr|-k)  +  -  JV  =  0, 


=  (») 

where  F(x)  =  Crrlh  —  t^rA-h  is  the  difference  between  external  and  internal 
normal  tractions  acting  on  the  shell. 

Integrating  (10)  from  —h  to  h  in  the  z  direction  and  using  (1)  and  (7) 
we  have: 

dN 

RiOrtU  -  +  -R-a-  +  Q  =  0, 
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where  T{x)  =  <rr6\h  —  tfre\-k  is  the  difference  between  external  and  internal 
shear  tractions  acting  on  the  shell. 

Multiplying  (10)  by  z  and  integrating  from  —h  to  h  in  the  z  direction 
and  using  (1)  and  (7)  we  have: 

/■*  ,  dcre  ,  daee  ,  „  X  .  dM  tn  ,  \ 

The  last  integral  in  this  equation  is: 

{R->rz)zarB\-h-  f  iR+2z)aredz  ts  -RQ  +  (R+h)harB\h  +  {R-h)hare\-k- 

J 

Using  these  formulas  we  arrive  at  the  equation  that  includes  the  bending 
moment  M: 

^  -  e  =  -5(1).  (13) 

where  S{x)  =  (l  +  h/.R)/iar«U  +  (l  -h/i2)h<rre|-fc  is  the  ratio  of  the  moment 
due  to  external  and  internal  shear  tractions  to  the  local  radius  of  curvature. 

Substituting  (8)  into  equations  (11)-(13)  we  arrive  at  the  following  equa¬ 
tions  in  terms  of  displacements: 


i?  h  Eh  ' 


^  +  “'T'T<-r 

«!-./)_  r(i-r=) 

2R  2Eh  ’ 


h(f>xx 


3(1  -  v) 


Zhx<l>x  +  ti' 


,3(1 -t/)  .3(1-*/)  3  5(1-*/^) 


*  2h  '  '  2Rh  ^  2h  2  Eh^  ■ 

Subscripts  “x"  and  “ii”  in  these  equations  indicate  partial  derivatives  from 
functions  u®,  w°  and  <l>. 

We  introduce  normalized  variables  u,  w  and  (  as: 


w°  =  Lw]  v°  =  Xu;  X  =  Xf. 
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Substituting  these  formulas  into  (14)  we  have: 

3  —  1/ 

ix;«  +  -  «€«(0  +  “«(0K0  +  MO  =  -/(Oi 


««  +  a(0«f  +  — ^«+  (15) 


=  -<({); 

+  3a(0^<  - 


where 


d(Infe)  d  R  L  h{0 

<0  ,K0  '*(«)  I  ' 


_  nO(i  +  t/)z  _  r(0(i- 1/^)1  _  „ 

J\0  ru/ts  y^\0  ~  nT<u/c\ 


3  5(0(1-1/^)!^ 
2  £;h3(0  • 


i:/i(£)  2EhiO 

This  system  of  -ordinary  differentiaJ  equations  of  the  second  order  can  be 
solved  with  the  foUowing  boundary  conditions: 


wiO)  =  u(0)  =  ^(0)  =  wil)  =  «(1)  =  <^(1)  =  0.  (16) 

Conditions  (16)  express  the  requirement  of  continuity  of  displacements  along 
a  shell  of  a  closed  cross-section.  In  addition,  we  fix  the  displacements  at  the 
origin  of  the  x  coordinate  system  setting  them  equal  to  zero.  This  is  a  usual 
condition  of  fixing  a  point  in  an  elastic  body  in  order  to  predefine  its  rigid 
translation. 

Equations  (15)  with  boundary  conditions  (16)  allow  us  to  find  the  dis¬ 
placement  of  the  shell  in  the  x  coordinate  system.  The  x  coordinate  is 
counted  along  the  cross-section  of  the  middle  surface  of  the  shell;  the  direc¬ 
tion  of  the  displacement  depends  on  the  local  geometry  of  the  middle  surface 
as  tn  is  perpendicular  to  the  middle  surface  and  u  is  tangent  to  it. 


The  Shape  of  the  Shell 

In  order  to  relate  the  displacements  of  the  shell  in  the  x  coordinate  system 
to  its  geometrical  configuration,  we  introduce  a  reference  rectangular  coor¬ 
dinate  system  (zq,  yo)  in  the  plane  of  deformation  (Fig.  3).  The  geometry  of 
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the  middle  surface  of  the  shell  can  be  given  in  the  (xq,  j/o)  coordinate  system 
in  a  parametric  form: 


xo  =  xo(t),  yo  =  yo(’'), 


(17) 


where  parameter  r  spans  an  interval  0  <  r  <  T.  This  parameter  can  be 
chosen,  for  example,  as  a  central  angle  a  changing  between  0  and  2ir. 

The  length  £  and  the  local  curvature  radius  Jt(r)  can  be  found  as: 


i  =  r  +  in  R(t)  =  (18) 

Jo  xoVo  —  yoxo 


xoVo  —  VoXo 

The  local  coordinate  x  is  related  to  the  parametric  equations  (17)  as: 


=  \/ifT7gdr. 


(19) 


Normalized  coordinate  xfL  is  calculated  using  (18)  and  (19), 

The  process  of  calculating  the  deformation  of  a  shell  under  a  given  load¬ 
ing  includes  the  numerical  solution  of  equations  (15)  with  boundary  condi¬ 
tions  (16),  and  consequent  calculation  of  the  shape  of  the  shell  after  defor¬ 
mation  in  a  reference  coordinate  system  (loij/o)*  The  latter  procedure  has 
simple  numerical  implementation  based  on  the  fact  that  displacement  u  is 
tangent  to  the  middle  surface  of  the  shell  and  w  is  normal  to  this  surface. 


Numerical  Examples 

In  this  section  we  give  two  numerical  examples  considering:  (1)  the  de¬ 
formation  of  a  circular  hollow  grain  with  and  without  cracks;  and  (2)  the 
deformation  of  an  elliptical  hollow  grain.  In  both  cases  grains  deform  due 
to  compression  along  the  vertical  axis.  The  magnitudes  of  acting  forces  are 
highly  exaggerated  in  order  to  make  the  deformation  visible  in  the  plots  pre¬ 
sented  below.  In  reality,  the  forces  of  such  magnitudes  will  result  in  plastic 
yield  of  grains. 


Circular  Grain 

We  examine  the  deformation  of  a  circular  hollow  grain  with  and  without 
radial  cracks  (Fig.  4)  due  to  compression  along  the  vertical  direction.  Two 
radial  fractures  dramatically  reduce  the  stiffness  of  the  grain  (Fig.  5).  This 
effect  may  be  used  to  explain  the  existence  of  the  shallow  low  velocity  zones 
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in  pillow  basalts.  The  closure  of  cracks  due  to  confining  pressure  increasing 
with  depth  will  result  in  a  sharp  stiffness,  and  seisnaic  velocity  increase. 

Normal  (w)  and  tangential  (u)  displacements  along  the  surface  of  the 
circular  hollow  grains  are  plotted  in  Fig.  6.  Displacements  in  the  fractured 
grain  (curves  a  and  c)  are  significantly  larger  than  in  the  grain  without 
cracks  (curves  b  and  d). 

Elliptical  Grain 

Computations  show  that  the  stiffness  of  an  unfractured  elliptical  hollow 
grain  (Fig.  7)  in  the  vertical  direction  is  close  to  the  stiffness  of  an  un- 
fractured  circular  grain.  However  the  stiffness  of  the  elliptical  grain  in  the 
horizontal  direction  is  more  than  twice  its  vertical  stiffness.  The  deformed 
shape  of  the  vertically  loaded  elliptical  grain  (a)  is  compared  to  its  initial 
shape  (b)  in  Fig.  8. 

Conclusions 

We  presented  a  theory  of  the  two-dimensional  deformation  of  hollow  grains 
under  external  loading,  modeling  them  as  cylindrical  shells  of  a  closed  cross- 
section.  This  theory  allows  us  to  describe  the  deformation  of  arbitrarily- 
shaped  grains  with  walls  of  varying  thickness,  and  external  and  internal 
fractures.  Fractures  dramatically  reduce  the  stiffness  of  the  grains.  This 
fact  may  be  used  to  explain  the  existence  of  the  shallow  low  velocity  zones 
in  pillow  basalts.  The  closure  of  cracks  due  to  confining  pressure  increas¬ 
ing  with  depth  will  result  in  a  sharp  stiffness,  and  seismic  velocity  increase. 
The  theory  presented  in  this  paper  can  be  directly  applied  to  estimating 
the  acoustical  and  mechanical  properties  of  basaltic  piUows  and  pelagic  sed¬ 
iments. 
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Figure  1:  A.  Arbitrarily-shaped  hollow  cylindrical  shell  and  a  local  rect¬ 
angular  coordinate  system  B.  local  coordinate  systems  {€,r)  and 


Figure  2:  A.  The  assumed  pattern  of  the  shell’s  deformation.  B.  forces  and 
moments  in  the  shell. 


Figure  3:  The  shell  in  a  reference  rectangular  coordinate  system  (xo,  yo)- 


Figure  4:  Circular  hollow  grain  with  two  radial  cracks. 


Figure  6;  The  deformed  shapes  of  circular  hollow  grains  without  cracks  (a) 
and  with  two  cracks  (b)  due  to  vertical  compression.  Curve  c  represents 
the  initial  circular  shape. 


Figure  6:  Displacements  along  the  surface  of  the  circular  hollow  grains: 
displacement  w  in  the  fractured  (a)  and  unfractured  (b)  grain;  displacement 
u  in  the  fractured  (c)  and  unfractured  (d)  grain. 


'  •  ‘  '1^^  111  ■  ^**"^1  I  I  I 
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Figure  8:  The  deformed  shape  of  the  elliptical  hollow  grain  (a)  compared 
to  its  initial  shape  (b). 


Pillow  Basalts  Compaction  due  to  Grain 

Collapse 


Abstract 

A  possible  mode  of  compaction  of  pillow  basalts  and  pelagic  sediments  (both 
formed  by  hollow  grains)  is  by  the  collapsing  of  the  grains  due  to  increas¬ 
ing  confining  pressure.  We  derive  simple  criteria  for  the  collapse  of  2-D 
thick-walled,  and  thin-wahed  circular  grains  subject  to  external  and  inter¬ 
nal  pressure. 

Introduction 

Assuming  that  the  compaction  of  pillow  basalts  and  pelagic  sediments  may 
occur  as  the  collapse  of  hollow  grains  that  form  those  rocks,  we  derive  simple 
criteria  for  the  compaction  due  to  increasing  confining  pressure. 

We  examine  the  deformation  of  a  2-D  thick-walled  cylinder  due  to  in¬ 
ternal  and  external  pressure.  To  find  the  criterion  of  the  collapse  of  such 
a  cylinder  we  employ  the  solution  for  the  elasto-plastic  deformation  of  a 
thick-walled  tube.  The  collapse  occurs  if  the  entire  cylinder  is  in  the  mode 
of  plastic  deformation.  In  this  solution  we  use  the  Mises  plastic  criterion. 

In  addition  we  give  simple  formulas  for  the  collapse  of  a  thin-walled 
cylindrical  tube  (shell). 

These  formulas  can  be  used  for  estimating  the  depth  of  transition  from 
low-velocity  to  high-velocity  zones  in  the  ocean  bottom  rocks.  This  transi¬ 
tion  is  presumably  connected  with  the  compaction  of  hollow  grains. 

Deformation  and  Collapse  of  a  Thick- Walled  Shell 

In  this  section  we  follow  the  solution  for  the  elasto-plastic  deformation  of  a 
thick-walled  tube  (e.g.  Godfrey,  1969). 

The  equation  of  bcilance  of  a  thick-walled  tube  subjected  to  the  internal 
pressure  Pi  and  external  pressure  P2  (Fig.  1)  in  the  cylindrical  coordinate 
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system  {r,0)  is: 


da  TT  =  (1) 

dr  r 

where  agg  and  ffrr  are  normal  stresses.  Equation  (1)  holds  both  for  plastic 
and  for  elastic  zones  of  the  tube. 

The  solution  of  (1)  with  the  boundary  conditions  Crr  =  Pi  at  r  =  a  and 
Crr  =  P2  at  r  =  h  for  the  elastic  deformation  of  the  tube  is: 


agg  = 


Pia^-Pib^  a^{^(P2-Pi)l 


52 -a2 


P^a'^-P^b'^  ,  aH\P2-Pi)l 

~  b^-a^  62  _  q2 

where  a  is  the  internal  and  b  is  the  external  radius  of  the  tube. 

The  Mises  plasticity  condition  in  the  cylindrical  coordinate  system  is: 


F  =  krr  - 


=  2*, 


where  k  is  the  plastic  limit  under  pure  shear  conditions,  and  k  =  oo/>/3, 
where  <to  is  the  plastic  limit  under  pure  tension. 

Substituting  (2)  and  (3)  into  (4)  we  find  that  if  P2  >  Pi 

^  2aH^{P2-Pi)l 
b^-a^  r2* 

The  function  P  increases  with  decreasing  r.  This  means  that  the  plasticity 
zone  will  be  initiated  at  the  internal  surface  of  the  tube  if  p  =  2k,  or 


P2  -  Pi  = 


k(b^  -  a2) 


As  the  difference  between  P2  and  Pi  increases,  the  plastic  zone  propa¬ 
gates  through  the  tube  towards  its  outer  surface.  The  whole  tube  will  be  in 
the  plastic  state  if  the  plastic  zone  reaches  the  external  radius  b. 
Substituting  (4)  into  (1)  we  find  that  in  the  plastic  zone 


Solving  (5)  with  the  boundary  condition  <7rr  =  —Pi  at  r  =  a  we  find 


Crr  =  -2/:  In - Pi. 

a 
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The  plastic  zone  reaches  r  =  6  if  Orr  =  —P2  r  =  b.  Substituting  this 
condition  into  (6)  we  find  that  this  happens  if 

P2-P1  =  2Jbln-.  (7) 

a 

Criterion  (7)  can  be  used  for  estimating  the  confining  pressure  causing 
the  collapse  of  a  hollow  grain. 

Given  k=  2  •  10®  Pa  and  6/a=l.l,  we  find  that  collapse  will  occur  at  the 
differential  pressure  P2  —  Pi  =380  bar. 

Collapse  of  a  Thin* Walled  Shell 

Normal  stress  acting  parallel  to  the  surfaces  of  a  thin-walled  cylindrical  shell 
under  differential  pressure  P  can  be  easily  found  as: 


PR 


where  R  is  the  radius  and  h  is  the  thickness  of  the  shell  (see  “The  Mechanics 
of  Hollow  Grains”). 

Thus  the  shell  collapses  if 

p  =  =  Vikj.  (8) 

Criterion  (8)  is  slightly  different  from  (7).  The  latter  has  the  following 
approximate  expression  as  a  approaches  b: 

P  =  2fcln(6/a)  w  2k{b—  a)/ a  —  2khla. 

This  difference  is  due  to  the  thin  shells  theory  approximations. 

Conclusions 

Assuming  that  the  compaction  of  pillow  basalts  and  pelagic  sediments  (both 
formed  by  hollow  grains)  is  the  result  of  the  collapse  of  the  grains  due  to 
increasing  confining  pressure,  we  derived  simple  criteria  for  the  collapse  of 
2-D  thick-walled,  and  thin- walled  circtilar  grains.  These  criteria  can  be  used 
for  estimating  the  depth  of  transition  between  the  low-velocity  and  the  high 
velocity  zones. 

Reference 
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Elasto-plastic  deformation  of  a  thick- walled  shell. 


Modeling  the  Seismic  Properties  of  Shallow  Oceanic 
Crust:  Effects  of  Extrusive  Morphology  on  Seismic 

Velocity 

D.  Moos.  J.  Dvorkin,  (both  at:  Dept,  of  Geophysics,  Stanford  U., 

Stanford,  CA  94305:  415  723-3464) 

Recent  seismic  data  indicate  that  seismic  velocities  at  the  top  of  very 
young  Oceanic  Layer  2  basalts  are  Vp^.5  km/s  and  Vs  =  0.6-0.8 
km/s.  The  low  velocity  zones  (LVZ's)  are  less  than  200  m  thick  and 
have  velocities  that  are  either  independent  of  depth  or  increase  very 
slowly.  Some  results  suggest  that  these  zones  are  absent  at  ridge 
crests,  increasing  in  thickness  with  crustal  age.  We  present  here  a 
series  of  micromechanical  models  to  explain  these  results  in  terms  of 
the  observed  characteristics  of  the  dominant  morphologies  of 
extrusive  basalts.  The  principle  result  of  these  exercises  is  that  there 
is  no  simple  relationship  between  seismic  velocities  and  porosity. 

These  models  rely  on  the  concept  of  pillow  basalts  and  basalt 
breccias  as  granular  aggregates  of  unit  elements.  The  velocity  of 
these  aggregates  is  dependent  on  the  stiffness  of  the  elements  and  of 
their  contacts.  This  approach  contrasts  with  effective  medium 
theories,  which  explain  velocity  variations  by  adding  to  a 
homogeneous  solid  material  a  set  of  isolated  cracks  with  a  well- 
defin^  porosity  and  distribution  of  aspect  ratios.  Early  results  to 
calculate  velocity  variations  by  varying  the  radius  of  contact  of 
pillows  modeled  as  solid  sphericd  elements  provided  a  good  fit  to  the 
velocity  data.  We  present  here  an  extension  of  the  m^el  in  which 
pillow  tubes  are  approximated  as  thin-walled  cylinders.  The  model 
can  accommodate  arbitrary  pillow  cross  sections,  and  can  include 
radial  (cooling)  cracks.  Depending  on  wall  thickness  and  the  number 
of  cracks,  pillow  tubes  modelled  in  this  way  can  have  very  small 
stiffnesses,  resulting  in  seismic  velocities  for  aggregates  composed 
of  these  elements  that  are  quite  small.  An  increase  in  contact  pressure 
accompanying  burial  has  very  little  effect  on  the  properties  of  these 
tubes;  however,  when  the  pressure  exceeds  the  strength  of  the 
cylinders,  they  will  collapse.  At  depths  corresponding  to  the  collapse 
pressure,  a  rapid  increase  in  velocity  may  result.  The  results  agree 
qualitatively  with  seismic  data.  Coupled  with  earlier  attempts  to 
model  contact  stiffnesses  of  pillows,  these  results  extend  our  efforts 
to  understand  the  velocity  structure  of  the  shallow  oceanic  crust  in 
terms  of  the  properties  of  the  dominant  morphological  types. 
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Abstract 

Prediction  of  wave  propagation  in  a  submarine  environment  requires  mod¬ 
elling  the  acoustic  response  of  ocean  bottom  sediments  which  generally  con¬ 
sists  of  porous  granular  materials  partially  or  wholly  saturated  with  water. 
The  effect  of  anisotropy  has  to  be  incorporated  into  the  model  in  order  to 
simulate  more  realistic  responses. 

Following  Biot’s  theory,  we  present  a  formulation  for  seismic  wave  ve¬ 
locities  in  a  general  anisotropic  poroelastic  medium.  We  also  identify  the 
anisotropic  parameters  that  need  to  be  evaluated  or  measured  for  the  pur¬ 
pose  of  estimating  the  velocities. 

Introduction 

Poroelasticity  or  the  mechanics  of  porous  media  is  of  fundamental  impor¬ 
tance  in  a  variety  of  fields  and  has  received  a  great  deal  of  attention  in 
recent  times.  It  is  believed  to  play  a  major  role  in  rock  failure  phenom¬ 
ena  and  triggerring  of  natural  and  induced  seismicity;  in  proper  selection  of 
underground  waste  disposal  sites;  enhanced  oil  recovery;  and  in  underwa¬ 
ter  acoustics  involving  wave  propagation  in  water  saturated  porous  marine 
sediments  of  the  ocean  bottom. 

Fluid  saturation  affects  the  acoustic  properties  of  porous  media  by  tend¬ 
ing  to  stiffen  them  in  compression  and  to  either  stiffen  or  soften  them  in 
shear.  The  effects  depend  in  complicated  ways  on  the  geometry  and  stiff¬ 
ness  of  the  pore  space,  the  stiffness  and  viscosity  of  the  fluid,  the  density  of 
all  components  and  the  frequency  of  wave  propagation. 

The  more  complex  the  internal  structure  and  constituents  of  the  rocks, 
the  greater  is  the  number  of  parameters  needed  to  characterize  its  proper¬ 
ties.  For  the  most  simple  case  of  isotropic,  elastic,  non-porous  solid  only  two 
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independent  parameters  are  needed;  four  parameters  are  necessary  to  char¬ 
acterize  an  isotropic  porous  medium  while  for  anisotropic  porous  media  the 
number  increases  depending  upon  the  symmetry  class  of  the  anisotropy  with 
28  independent  parameters  being  required  for  the  most  general  anisotropic 
poroelastic  medium. 

Much  of  our  understanding  of  saturation  effects  and  poroelastic  phe¬ 
nomena  has  been  based  on  the  formulations  of  Biot  (1-6)  whose  equations 
have  formed  the  basis  for  solving  particular  problems  in  poroelasticity  and 
has  long  been  regarded  as  a  sort  of  standard.  A  great  deal  of  the  literature 
pertinent  to  the  study  of  sediment  acoustics  (e.g.  7-10)  has  however  focused 
on  the  isotropic  aspects  of  Biot's  theory  and  there  has  been  a  general  lack  of 
using  the  corresponding  anisotropic  theory  for  realistic  granular  sediments 
and  rocks.  We  have  begun  to  extend  these  descriptions  to  anisotropic  porous 
solids  starting  with  the  anisotropic  constitutive  equations  presented  by  Biot. 

Biot’s  Theory 

Biot’s  theory  is  a  promising  approach  for  modelling  acoustic  wave  propa¬ 
gation  in  ocean  sediments  which  generally  consist  of  elastic  or  viscoelastic 
porous  granular  materials  saturated  with  water.  It  is  a  phenomenological 
theory  with  no  restrictions  being  made  a  priori  regarding  the  geometry  and 
shape  of  cracks  and  pores.  The  fluid  constituent  is  modelled  as  an  inter¬ 
connected  compressible  continuum  within  the  porous  solid  matrix,  each  of 
them  capable  of  having  independent  motions.  Biot  considered  both  low  and 
high  frequency  limits  in  his  analysis  of  propagation  of  harmonic  waves  in 
an  unbounded  fluid-saturated  porous  elastic  medium  and  predicted  the  ex¬ 
istence  of  three  kinds  of  body  waves:  two  compressional  or  dilational  and 
one  rotational.  The  two  distinct  types  of  compressional  waves  are  called  the 
fast  and  s/ou'  waves  or  P-waves  of  the  first  and  second  kind. 

The  P-waves  of  the  first  kind  are  similar  to  those  observed  in  ordinary 
elastic  media  and  is  characterized  by  high  phase  velocity,  small  dispersion 
and  low  attenuation  due  to  the  fact  that  the  motions  of  the  porous  frame  and 
the  fluid  are  nearly  in  phase  and  hence  viscous  losses  are  relatively  small.  On 
the  other  hand,  compressional  waves  of  the  second  kind  are  diffusive  waves, 
have  low  phase  velocities,  large  dispersion  and  are  highly  attenuated  with 
the  frame  and  fluid  components  moving  largely  out  of  phase.  In  general  it  is 
quite  difficult  to  observe  the  slow  wave  in  real  situations,  but  the  cumulative 
effect  of  energy  conversion  from  fast  to  slow  waves  at  interfaces  between 
dissimilar  porous  media  and  between  a  fluid  and  a  porous  medium  (e.g. 
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ocean  bottom)  is  quite  signiAcant  and  leads  to  considerable  dissipation  of 
the  fast  P-wave. 

Anisotropic  Poroelasticity;  Derivations 

In  this  section  we  present  a  formulation  for  estimating  seismic  velocities  in 
general  anisotropic  porous  solids.  Following  Biot’s  theory  (5),  the  constitu¬ 
tive  equations  for  poroelasticity  in  the  most  general  case  of  anisotropy  may 
be  written  as  (summation  convention  implied): 

+  Af.jC  (1) 

=  (2) 

where  is  the  stress  tensor,  is  the  strain  tensor,  pj  the  fluid  pressure, 
Cijku  Mij  and  M  are  material  constants,  (  represents  the  increment  of  fluid 
content  and  is  given  by 


C=-V.w  (3) 

where 

w  =  ^(U  -  u)  (4) 

Here  XJ  is  the  average  fluid  displacement  vector,  u  the  displacement  of  the 
solid  matrix  and  ^  the  porosity. 

The  equations  of  motions  governing  wave  propagation  are  (5): 


Ciikt€ki,i  +  =  PUi  +  PjWi  (5) 

-  -  M C,i  =  4  mijWj  4  bijWj  (6) 

or 

C,jkivi,kj  -  MijWk,k}  =  4  pjWi  (7) 

-  4  Mwk,ki  =  Pj<ii  4  lij(wj)  (8) 


where  p  =  (1  -  4>)p,  •¥  <i>pj  ,  p,  and  pj  being  the  solid  and  fluid  density 
respectively.  m,j  and  6,j  are  the  anisotropic  virtual  mass  and  drag  coeffi¬ 
cients,  together  forming  the  viscodynamic  operator  X.,y.  These  six  equations 
for  the  unknown  vector  components  u,  and  u»,  govern  the  propagation  of 
waves  in  general  anisotropic  porous  media. 

First  consider  the  high  frequency  limit  which  is  the  same  as  the  case  of 
an  inviscid  fluid  since  for  large  frequencies  the  phase  velocity  with  viscous 
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diBsipation  tends  towards  that  foi  inviscid  fluid.  In  this  limit  the  term 
in  the  first  derivative  of  tr,  may  be  dropped  from  the  second  of  the  two 
equations  of  motions.  Take  plane  wave  solutions  of  the  form  u,-  =  j4,'/(r)  and 
Wi  =  Bif{T)  where  t  =  t  —  n,-  are  the  components  of  the  unit  vector 

along  the  direction  of  propagation  and  V  the  phase  velocity.  Substituting 
this  into  equations  (7)  and  (8)  after  dropping  the  viscous  dissipation  term 


we  get; 

CijkinjnkAi  -  MijTiknjBk  =  pV^Ai  +  pjV^Bi  (9) 

-  MijUiTiiAj  +  MrikniBk  -  PjV^Ai  +  mijV^Bj  (10) 

These  formulas  can  be  written  as: 

TiiA,  -  'fikBk  =  pV^Ai  +  pjV^B,  (11) 

-  'IfjiAj  +  6kiBk  =  p^y^Ai  +  mijV^Bj  (12) 

where  r,i  =  Cijkinjrik  is  the  Christoffel  tensor  (e.g.  11,  12),  7^*  =  MijTijTik 
and  6k\  —  A/ufcnj.  Solving  equation  (12)  for  B  gives; 

B  =  (7ii  +  P/V2)D-^A  (13) 

with  Du  =  6n  -  muV^.  Substituting  this  in  equation  (11)  we  get 

TiiAi  -  (7i;  +  PjV^)Heii  -  m<,V2)-U,  =  pV^Ai  (14) 

This  can  be  written  in  the  form  of  Christoffers  equation  as 

(fl/  -  pV%i)  Ai  =  0  (15) 

where  F,;  is  the  Poroelastic  Christoffel  Tensor  given  by: 


r.7  =  Ti,  -  (7ii  +  P/V^fieu  -  (16) 

The  roots  of  the  determinantal  equation 

|f7;  -  pVHn\  =  0  (17) 

give  the  velocities  (or  eigenvalues)  corresponding  to  the  A/ -eigen vectors. 
Now  if  we  solve  equation  (12)  for  A  we  get 


A,= 


(g.>  -  m,kV^) 
(7.,  +  PjV^) 


(18) 
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which  on  substituting  back  in  equation  (11)  gives 

{iik-PjVHik)Bk  =  0 

(19) 

where 

^ik  “  ,  T/2  f^k 

7ii  +  P/V^ 

Here  again  the  determinantal  equation 

(20) 

|7ifc  -  pfy%k\  =  0 

(21) 

gives  the  eigenvalues  corresponding  to  the  Hjt-eigenvectors. 

Considering  two  limiting  cases:  one  when  C  =  —divw  = 
of  motion  reduce  to: 

0,  the  equations 

Cijkmjcj  =  pui 

(22) 

(23) 

Again  substituting  plane  wave  solutions  of  the  form  m  =  Ai/(t— n,ii/V) 
we  get; 

iTi,-pV^6i,)A,:^0  (24) 

liiAi  =  -p/V^A,  (25) 

Solving  these  gives  an  Initial  value  for  which  can  be  used  in  equations 
(16)  and  (17). 

Another  limiting  case  is  when  pj  —  0.  Now  we  get 


ri,A,-^ikBk  =  pV^Ai-i-pfV^S, 
pjAi  +  mijSj  —  0 

These  two  equations  ultimately  give 


(26) 

(27) 


(f;;-pV2^i,)A,  =  0  (28) 

with 

f7,  =  Ti,  +  (29) 

mu 

So  far  the  derivation  was  without  fluid  viscosity  which  is  appropriate 
for  the  high  frequency  limit.  Now  talking  into  aiccount  fluid  viscosity  rj 
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aod  agdn  substituting  plane  wave  solutions  in  the  equations  of  motions  we 
get  a  frequency- velocity  (w  —  V)  relationship  expressed  in  the  form  of  a 
determinant  as: 

PfJV^iij  -f  Lf^ij  uV^rriij  —  —  u>6ij 

(30) 

Each  term  is  a  3  x  3  determinant  making  the  whole  a  6  x  6  one.  In  the 
above,  the  solid  skeleton  is  considered  to  be  perfectly  elastic  and  dissipation 
is  only  due  to  fluid  viscosity. 

Biot’s  Parameters 

The  basic  parameters  of  Biot’s  theory  can  be  divided  into  two  categories:  the 
passive  constants  (  Ciju,  M  )  which  depend  on  the  stiffnesses  and  bulk 
modulus  of  the  solid  and  fluid  components,  their  densities  and  the  porosity 
and  a  second  category  (  Lij  -  the  viscodynamic  operator  )  which  includes  a 
number  of  interrelated  parameters  dependent  on  the  relative  motion  of  the 
fluid  with  respect  to  the  skeletal  frame.  The  viscodynamic  operator  involves 
the  first  and  second  derivative  with  respect  to  time  of  the  variable  tn,-. 

Lij  =s  (^0^^  + 

The  coefficients  m,j  and  bij  depend  on  the  fabric  or  microstructure  of 
the  constituents  and  the  wave  frequency. 

The  fundamental  problems  associated  with  the  application  of  these  equa¬ 
tions  to  wave  propagation  consists  of  determination  of  the  various  moduli 
and  coefficients  in  terms  of  the  solid,  fluid  and  porous  media  properties  and 
comparision  with  suitable  experimental  data.  Perhaps  a  better  understand¬ 
ing  of  the  parameters  may  be  obtained  by  considering  wave  propagation  in 
isotropic  porous  media.  Then  the  equations  of  motion  simplify  to  (5): 


-KH  -  p)ejj,i  -  C C.i  =  pUi  •¥  Pj th,  (32) 

-  M(l,i  =  pjUi  tnwi  +  6u),  (33) 

It  has  been  shown  (13,  14)  that  the  coefficients  H,  C,  M  and  p  can 
be  expressed  in  terms  of  the  porosity,  the  properties  of  the  fluid  and  solid 
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constituents  and  the  elastic  moduli  of  the  dry  porous  medium.  The  results 
can  be  written  as  (15): 

M  =  K,/[l  -  Kb/K,  +  <f>iKJKj  -  1)]  (34) 

C  =  (1  -  Kb/K,)M  (35) 

E  =  il-  Kb/K,)C  +  iTi  +  (4/3)/x  (36) 

where  A'/  is  the  bulk  modulus  of  the  fuld,  Kt  that  of  the  solid  and  Kt  and 
are  the  bulk  and  shear  moduli  of  the  dry  porous  medium,  also  known  as 
the  skeletal  frame  moduli. 

In  a  few  cases  these  coefficients  or  some  equivalent  parameter  have  been 
determined  for  porous  material  by  carrying  out  jacketed  and  unjacketed 
hydrostatic  compression  tests  and  shear  tests  (16,  17,  18).  It  has  been 
suggested  (19),  that  the  moduli  of  the  skeletal  frame  Kt  and  fx  shoiild  be 
assumed  to  be  complex  to  account  for  viscous  effects. 

Apart  from  these  constitutive  constants  the  equations  of  motion  for  a 
fluid  saturated  permeable  solid  also  has  the  frequency  dependent  viscody- 
namic  operator  L(wi)  =  mwi  +  btiii  which  contains  the  inertial  drag  and 
viscous  dissipative  effects  due  to  the  pore  fluid  motion.  The  simplest  form 
for  m  and  6  as  originally  proposed  by  Biot  is  6  =  v/k  and  m  s=  Qpj/<f>  where 
V  is  the  dynamic  fluid  viscosity,  k  the  permeability  for  steady  state  flow  and 
a  is  a  structure  factor  for  the  pore  space.  This  simple  form  of  T  is  implicitly 
a  low  frequency  approximation  as  it  uses  the  steady  state  Darcy  term  plus  a 
first  order  inertial  correction  term.  Biot  presented  a  method  for  evaluating 
the  drag  coefficient  by  determining  the  average  velocity  of  liquid  in  a  circu¬ 
lar  cylinder  subject  to  an  oscillating  pressure.  A  closely  related  method  for 
evaluating  the  drag  and  virtual  mass  coefficients  for  specific  pore  geometries 
has  been  presented  (20)  wherein  the  coefficients  are  determined  by  solution 
of  a  boundary  value  problem  for  a  viscous,  compressible  fluid  in  a  single  pore 
of  specified  shape.  Using  finite  difference  or  finite  element  solutions  for  fluid 
displacement  it  is  now  feasible  to  determine  (for  the  isotropic  case)  the  drag 
and  virtual  mass  coefficients  in  Biot’s  equations  as  a  function  of  frequency 
for  more  complicated  pore  gf  imetries  (21). 
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Conclusions 

Biot’s  phesomenologicaJ  theory  for  wave  propagation  in  fluid  saturated  porous 
media  presents  a  promising  framework  for  modelling  the  interaction  between 
seismic  waves  impin^ng  on  the  sea  floor  and  the  submarine  sediments.  Much 
work  has  been  done  for  isotropic  porous  media.  However  to  simulate  more 
realistic  response  the  effect  of  anisotropy  has  to  be  taken  into  account. 

We  have  presented  a  mathematical  formulation  based  on  Biot’s  anisotropic 
constitutive  equations  for  estimating  seismic  velocities  in  a  general  anisotropic 
poroelastic  medium.  Two  cases  have  been  considered:  that  of  an  inviscid 
fluid  (which  is  appropriate  for  the  high  frequency  limit)  and  that  of  a  viscous 
fluid,  which  is  valid  for  all  frequencies.  In  order  to  solve  the  determinantal 
equations  for  the  wave  velocities,  one  needs  to  know  the  relevant  anisotropic 
rock  parameters,  the  poroelastic  material  constants  and  the  anisotropic  vir¬ 
tual  mass  and  drag  coefficients  which  together  form  the  viscodynamic  oper¬ 
ator. 

W’hile  measurements  and  estimates  of  Biot  parameters  (both  the  mate¬ 
rial  constants  and  the  viscodynamic  operator)  for  isotropic  material  have 
been  made,  the  corresponding  anisotropic  Biot  coefficients  remain  to  be 
determined  and  an  important  obstacle  to  the  comparision  of  theory  with 
acoustic  measurements  in  anisotropic  porous  media  still  remains. 
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